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Abstract 

DDSCAT. 5a is a freely available software package which applies the "discrete dipole ap- 
proximation" (DDA) to calculate scattering and absorption of electromagnetic waves by targets 
with arbitrary geometries and complex refractive index. The DDA approximates the target by 
an array of polarizable points. DDSCAT. 5a requires that these polarizable points be located 
on a cubic lattice. DDSCAT. 5al0 allows accurate calculations of electromagnetic scattering 
from targets with "size parameters" 27ra/A < 15 provided the refractive index m is not large 
compared to unity {\m — 1| < 1). 

The DDSCAT package is written in Fortran and is highly portable. The program supports 
calculations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, rectangular 
solids, finite cylinders, hexagonal prisms, etc.). Target materials may be both inhomogeneous 
and anisotropic. It is straightforward for the user to "import" arbitrary target geometries 
into the code, and relatively straightforward to add new target generation capability to the 
package. DDSCAT automatically calculates total cross sections for absorption and scattering 
and selected elements of the Mueller scattering intensity matrix for specified orientation of the 
target relative to the incident wave, and for specified scattering directions. 

This User Guide explains how to use DDSCAT. 5al0 to carry out electromagnetic scatter- 
ing calculations. CPU and memory requirements are described. 
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1 Introduction 



DDSCAT.Sa is a Fortran software package to calculate scattering and absorption of electro- 
magnetic waves by targets with arbitrary geometries using the "discrete dipole approximation" 
(DDA). In this approximation the target is replaced by an array of point dipoles (or, more 
precisely, polarizable points); the electromagnetic scattering problem for an incident periodic 
wave interacting with this array of point dipoles is then solved essentially exactly. The DDA 
(sometimes referred to as the "coupled dipole approximation") was apparently first proposed 
by Purcell & Pennypacker (1973). DDA theory was reviewed and developed further by Draine 
(1988), Draine & Goodman (1993), and recently reviewed by Draine & Flatau (1994) and Draine 
(2000). 

DDSCAT.Sa is a Fortran implementation of the DDA developed by the authors. The pre- 
vious version, DDSCAT.5a9, was released 1998 December 15. The current version, DDSCAT.BalO, 
released 2000 June 15, adds a new "multisphere" target option. DDSCAT.Sa is intended to 
be a versatile tool, suitable for a wide variety of applications ranging from interstellar dust to 
atmospheric aerosols. As provided, DDSCAT.SalO should be usable for many applications 
without modification, but the program is written in a modular form, so that modifications, if 
required, should be fairly straightforward. 

The authors make this code openly available to others, in the hope that it will prove a useful 
tool. We ask only that: 

• If you publish results obtained using DDSCAT, please consider acknowledging the source 
of the code. 

• If you discover any errors in the code or documentation, please promptly communicate 
them to the authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public Li- 
cense) of the Free Software Foundation: you may copy, distribute, and/or modify the 
software identified as coming under this agreement. If you distribute copies of this soft- 
ware, you must give the recipients all the rights which you have. See the file doc/ copyleft 
distributed with the DDSCAT software. 

We also strongly encourage you to send email to the authors identifying yourself as a user of 
DDSCAT; this will enable the authors to notify you of any bugs, corrections, or improvements 
in DDSCAT. 

The current version, DDSCAT.SalO , uses the DDA formulae from Draine (1988). The 
code incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & Flatau 1991). 
The "lattice dispersion relation" (LDR) prescription (Draine & Goodman 1993) is used for 
determining dipole polarizabilities. 

This User Guide assumes that you have already obtained the Fortran source code for 



DDSCAT.SalO cither via the World Wide Web ( http: //www, astro .princeton. edu/ -^draine) 
or via anonymous ftp following the instructions in the README file.[J We refer you to the list of 
references at the end of this document for discussions of the theory and accuracy of the DDA 
[first see the recent reviews by Draine and Flatau (1994) and Draine (2000)]. In §| we describe 
the principal changes between DDSCAT.SalO and the previous releases.]^ The succeeding 
sections contain instructions for: 



^To obtain the README file: (1) anonymous ftp to astro.princeton.edu , (2) cd draine/scat/DDA/ver5a, and 
(3) get README. 

^ The first "official release" was DDSCAT. 4b, although DDSCAT. 4c - while never announced - was made 
available to a number of interested users. DDSCAT. 5a8 was released in 1997 April. DDSCAT. 5a9 was released 
in 1998 December. DDSCAT.SalO was released 2000 June 15. 
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• compiling and linking the code; 

• running a sample calculation; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations; 

• specifying target orientation; 

• changing the DIMENSIDNing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a UNIX system; 
slight changes will be necessary for non-UNIX sites, but they are quite minor and should present 
no difficulty. 

Finally, this User Guide can be obtained by any of the following means: 



• http : //xxx ■ lanl . gov/ abs/astro-ph/0008151 - you will be offered the options of down- 
loading 

— Latex source 

— Postscript 

— Other formats - click on this to obtain the UserGuide as a PDF file. 

• anonymous ftp from astro . princeton . edu, subdirectory draine/ scat/ddscat/ver5al0. 

— For postscript, get either the uncompressed postscript file UserGuide.ps or the 
gzipped file UserGuide. ps.gz (remember to set the binary option in ftp before 
using the get command). 

— Alternatively, you can get the LaTeX source UserGuide . tex . gz and the figures 
f l.eps.gz, f2.eps.gz, f3.eps.gz, f4.eps.gz, and f5.eps.gz. 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of 
the target, being limited only by the need to use an interdipole separation d small compared to 
(1) any structural lengths in the target, and (2) the wavelength A. Numerical studies (Draine 
& Goodman 1993; Draine & Flatau 1994; Draine 2000) indicate that the second criterion is 
adequately satisfied for calculations of total cross sections if 

\m\kd<l , (1) 

where m is the complex refractive index of the target material, and k = 27r/A, wh ere A is the 
wavelength in the surrounding medium (normally taken to be vacuum, but see §^]^). However, 
if accurate calculations of the scattering phase function (e.g., radar or lidar cross sections) are 
desired, a more conservative criterion 

|m|fcd<0.5 (2) 

will ensure that differential scattering cross sections dCgca/dfl are accurate to within a few 
percent of the average differential scattering cross section Csca/47r (see Draine 2000). 

Let V be the target volume. If the target is represented by an array of N dipoles, located 
on a cubic lattice with lattice spacing d, then 

V = Nd^ . (3) 



5 



We characterize the size of the target by the "effective radius" 

fleff = (3F/47r)i/3 , (4) 

the radius of an equal volume sphere. A given scattering problem is then characterized by the 
dimensionless "size parameter" 

X = fcocff = — - — . (5) 
The size parameter can be related to N and Imj/cd: 

1/3 



27raeff 62.04 / N 



A |m| VlO^ 

Equivalently, the target size can be written 

A /iV^^/^ 



I fed . (6) 



m 



fleff = 9.873— ( — ) • \m\kd . (7) 



Practical considerations of CPU speed and computer memory currently available on scientific 
workstations typically limit the number of dipoles employed to < 10^ (see §^for limitations 
on TV due to available RAM); for a given A^, the limitations on \m\kd translate into limitations 
on the ratio of target size to wavelength. 

For calculations of total cross sections Cabs and Csca, we require jml/cd < 1: 

A / N 62.04 / A^ ^ 



For scattering phase function calculations, we require Imj/cd < 0.5: 

A / N 31.02 / N ^^^^ 



m 



«e.<4.94-(^ or^<--ij^] . (9) 



It is therefore clear that the DDA is not suitable for very large values of the size parameter x, 
or very large values of the refractive index to. The primary utility of the DDA is for scattering by 
dielectric targets with sizes comparable to the wavelength. As discussed by Draine & Goodman 
(1993), Draine & Flatau (1994), and Draine (2000), total cross sections calculated with the 
DDA are accurate to a few percent provided A^ > lO'* dipoles are used, criterion (P is satisfied, 
and |to - 1| < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs. which show overall 
scattering and absorption efficiencies as a function of wavelength for different discrete dipole 
approximations to a sphere, with A^ ranging from 304 to 59728. The DDA calculations assumed 
radiation incident along the (1,1,1) direction in the "target frame". Figs. show the scat- 
tering properties calculated with the DDA for x ~ ka = 7. Additional examples can be found 
in Draine & Flatau (1994) and Draine (2000). 

3 DDSCAT.Sa 

3.1 What Does It Calculate? 

DDSCAT.Sa solves the problem of scattering and absorption by an array of polarizable point 
dipoles interacting with a monochromatic plane wave incident from infinity. DDSCAT.5a has 
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1 2 3 4 5 6 7 8 9 10 11 12 13 
x = ka 



Figure 1: Scattering and absorption for a sphere with m = 1.33 + O.Oli. The upper panel shows the 
exact values of Qsca and Qabs; obtained with Mie theory, as functions of x = ka. The middle and 
lower panels show fractional errors in Qsca and Qabsj obtained using DDSCAT with polarizabilities 
obtained from the Lattice Dispersion Relation, and labelled by the number of dipoles in each 
pseudosphere. After Fig. 1 of Draine & Flatau (1994). 




0123456789 10 11 
x = ka 



Figure 2: Same as Fig. ||, but for m = 2 + i After Fig. 2 of Draine & Flatau (1994). 
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Figure 3: Differential scattering cross section for m = 1.33 + O.Oli pseudosphere and ka = 7. Lower 
panel shows fractional error compared to exact Mie theory result. The N = 17904 pseudosphere 
has \m\kd = 0.57, and an rms fractional error in da/dO, of 2.4%. After Fig. 5 of Draine Sz Flatau 
(1994). 




Figure 4: Same as Fig. |3| but for m = 2 -\- i. The N 
an rms fractional error in da/dfl of 6.7%. After Fig. 



59728 pseudosphere has \m\kd 
of Draine & Flatau (1994). 



0.65, and 
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the capability of automatically generating dipole array representations for a variety of target 
geometries (see §19) and can also accept dipole array representations of targets supplied by the 
user (although the dipoles must be located on a cubic lattice) . The incident plane wave can have 
arbitrary elliptical polarization (see §^), and the target can be arbitrarily oriented relative to 
the incident radiation (see §[l^). The following quantities are calculated by DDSCAT.5a : 

• Absorption efficiency factor Qabs = Cabs/7''aeff' where Cabs is the absorption cross section; 

• Scattering efficiency factor Qsca = C'sca/'ra^g, where Csca is the scattering cross section; 

• Extinction efficiency factor Qcxt = Qsca + Qabs] 

• Phase lag efficiency factor Qpha, defined so that the phase-lag (in radians) of a plane wave 
after propagating a distance L is just ?T.tQpha7''a^ffi, where rit is the number density of 
targets. 

• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering prop- 
erties of the target for scattering directions specified by the user. 

• The radiation force efficiency vector Qp,. (see §p^). 

• The radiation torque efficiency vector Qr (see alq). 



3.2 Application to Targets in Dielectric Media 

Let to be the angular frequency of the incident radiation. For many applications, the tar- 
get is essentially in vacuo, in which case the dielectric function e which the user should sup- 
ply to DDSCAT is the actual complex dielectric function etargct(w), or complex refractive index 
^target {^) = Retarget of the target material. 

However, for many applications of interest (e.g., marine optics, or biological optics) the 
"target" body is embedded in a (nonabsorbing) dielectric medium, with (real) dielectric function 
emedium(w). Or (real) refractive index m^^cdiMmiy) = V^modium- DDSCAT. 5a is fully applicable 
to these scattering problems, except that: 

• The "dielectric function" or "refractive index" supplied to DDSCAT. 5a should be the 
effective dielectric function 

^target 



e(^) = --^-y- (10) 



'^medium V 



or effective refractive index: 



m{u:) = '^^-^^-^^';\ _ (11) 

'71medium(w) 



The wavelength A specified in ddscat .par should be the wavelength in the medium: 

A= ^""^ , (12) 

^medium 

where X^^ac — Zttc/cj is the wavelength in vacuo. 

The absorption, scattering, extinction, and phase lag efficiency factors Qabs, Qsca, and Qoxt cal- 
culated by DDSCAT will then be equal to the physical cross sections for absorption, scatter- 
ing, and extinction divided by Tra^fj - e.g., the attenuation coefficient for radiation propagating 
through a medium with a density rit of scatterers will be just a = ntQcxtT^alff. Similarly, the 
phase lag (in radians) after propagating a distance L will be utQphi 



2 

off- 
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The elements Sij of the 4x4 MueUer scattering matrix S calculated by DDSCAT will be 
correct for scattering in the medium: 



27rr 



S • Ii„, (13) 



where lin and Isca are the Stokes vectors for the incident and scattered light (in the medium), 
r is the distance from the target, and A is the wavelength in the medium (eq. |l^). See §^ for 
a detailed discussion of the Mueller scattering matrix. 

The time-averaged radiative force and torque (see § [l5| ) on a target in a dielectric medium 

are 

Frad = QprTTacff^rad , (14) 
Trad = Qr^raeff'^rad:^ , (15) 

where the time-averaged energy density is 

'^rad — ^medium 7^ , (1^) 

where Eq cos(cjt + cj)) is the electric field of the incident plane wave in the medium. 



4 What's New? 
4.1 DDSCAT.5a 

DDSCAT. 5a differs from previous versions in five major respects: 

1. Use of the new Generalized Prime Factor Algorithm (GPFA) developed by Clive Tem- 
perton (1992) for FFT calculations. The GPFA algorithm is generally faster than the 
previous algorithms, yet requires no more memory than the algorithm of Brenner (1969). 
See ^ 

2. Availability of a new algorithm for iterative solution of the system of complex linear 
equations. This algorithm is often faster than the algorithm of Petravic and Kuo-Petravic 
(1979) which was used through DDSCAT. 4b (and which remains available as an option 
in DDSCAT.Sa). See §|4[ 

3. Automatic calculation of the transverse components of the radiative force on the target. 



See ilb 



Capability to compute the electromagnetic torque on the target, due to absorption and 
scattering of light from the incident beam, as described by Draine and Weingartner (1996). 



See §15 



5. Automatic calculation of the 4x4 Mueller scattering matrix. See § |25| . 

6. Improved output handling, including FORTRAN unformatted binary write option as well 
as a net CDF interface and IDL postprocessing code. 

7. Finally, this much-improved and expanded User Guide! 

We also call users' attention to a minor but possibly confusing change: DDSCAT. 5a uses a 
different convention for specifying the target axes ai and 3.2 for certain target choices (RCTNGL, 
ELLIPS, ...). We think the new convention is more straightforward. See the discussion in § p^ 
below for details. 
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4.2 DSCAT.5a9 



DDSCAT.5a9 differs from DDSCAT.5a8 in one respect: 

• DDSCAT.5a8 contained an error in the calculation of the elements of the scattering 
matrix (see § p3| ) for scattering planes with 4> ^ Q. While DDSCAT.5a8 computed the 
scattering matrix correctly for (j)s = 0, for irregular targets only the 5ii clement was 
computed correctly for scattering planes with (ps 7^ 0. This bug has been corrected in 
DDSCAT.5a9. 



4.3 DDSCAT.SalO 

DDSCAT.5alO differs from DDSCAT.5a9 in two respects: 

• The public-domain LAPACK code has been replaced by the latest routines obtained from 



NETLIB ( jittp : / /www . net lib . org| ). Subroutine PRINAXIS did not function properly 
when the older LAPACK routines were compiled using the g77 compiler, most likely due 
to some variable not being SAVE-d. The new LAPACK code appears to execute properly. 

A new shape option - NSPHER - has been added, permitting a target to be defined as 
the union of the volumes of an arbitrary number of spheres. The locations and radii of 
the NSPH spheres are input via a file. Adding this option required the argument list of 
subroutine TARGET to be changed, and modifications were also made to subroutines REAPAR 
and REASHP. 

A new shape option - PRISMS - has been added (on 2002.02.12), enabling construction of 
a triangular prism target. Concomitant modifications were made to subroutines TARGET 
and REAPAR. 



5 Obtaining the Source Code 

The easiest way to download the source code is from the DDSCA T home page: 

http : //www ■ astro . princeton . edu/ ^draine/DDSCAT . html 
where you can obtain the file ddscatSalO . tar . gz - a "gzippcd tarfile" containing the com- 
plete source code and documentation. This can also be obtained by anonymous ftp from 
astro.princeton.edu, sudirectory draine/scat/ddscat/ver5al0. Note that it is a com- 
pressed (binary) file. 

The source code can be installed as follows. First, place the file ddscatSalO.tar .gz in the 
directory where you would like DDSCAT to reside. You should have at least 5 Mbytes of disk 
space available. 

If you are on a Linux system, you should be able to type 

tar xvzf ddscatSalO.tar .gz 
which will "gunzip" the tarfile and then "extract" the files from the tarfile. If your version of 
"tar" doesn't support the "z" option (e.g., you are running Solaris) then try 

zcat ddscatSalO .tar .gz | tar xvf - 
If neither of the above work on your system, try the two-stage procedure 

gunzip ddscatSalO . tar . gz 

tar xvf ddscatSalO . tar 
(The disadvantage of the two-stage procedure is that it uses more disk space, since after the 
second step you will have the uncompressed tarfile ddscatSalO .tar - about 3.8 Mbytes - in 
addition to all the files you have extracted from the tarfile - another 4.6 Mbytes). 
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Any of the above approaches should create subdirectories src, doc, misc, and IDL. The 
source code will be in subdirectory src, and documentation in subdirectory doc. 

If you are running Windows on a PC, you will need the "winzip" program, which can 
be downloaded from http : //www. winzip. com. winzip should be able to "unzip" the gzipped 
tarfile ddscatBalO . tar . gz and "extract" the various files from it automatically. 



6 Compiling and Linking 

In the discussion below, it will be assumed that the source code for DDSCAT.Sa has been 
installed in a directory DDA/ src . To compile the code on a Unix system, position yourself 
in the directory DDA/src. The Makefile supplied has compiler options appropriate for Sun 
Fortran under Solaris 2.x. If you have a different Fortran compiler, you will probably need to 
edit Makefile to provide the desired compiler options. Makefile contains samples of compiler 
options for selected operating systems, including HP AUX, IBM AIX, and SGI IRIX operating 
systems. 

DDSCAT can be compiled on Windows systems using standard Fortran compilers (e.g., 
Compaq or Microsoft Visual Fortran), and run successfully. Note, however, that many of the 
instructions below are specifically for Unix or Linux operating systems - some steps may have 
to be done differently on a non-Unix platform. 

So far as we know, there are only two operating-system-dependent aspects of DDSCAT. 5a: 
(1) the device number to use for "standard output", and (2) the TIMEIT routine. There is, 
in addition, one installation-dependent aspect to the code: the procedure for linking to the 



netCDF library (see §10.3 for discussion of the possibility of writing binary files using the 



machine-independent netCDF standard). 



6.1 Device Numbers IDVOUT and IDVERR 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error 
messages" , respectively. Normally these would both be set to the device number for "standard 
output" (e.g., writing to the screen if running interactively). The variables IDVOUT and IDVERR 
are set by DATA statements in the "main" program DDSCAT . f and in the output routine WRIMSG 
(file wrimsg.f ). Under Sun Fortran, DATA lDVDUT/0/ results in unbuffered output to "standard 
output" ; unbuffered output (if available) is nice so that if you choose to direct your output to a 
file (e.g., using ddscat >& ddscat.out &) the output file will contain up-to-date information. 
Other operating systems or compilers may not support this, and you may need to edit DDSCAT . f 
to change the two data statements to DATA lDVOUT/6/ and DATA lDVERR/6/, and edit wrimsg . f 
to change DATA lDVDUT/0/ to DATA lDVOUT/6/. 

6.2 Subroutine TIMEIT 

The only other operating system-dependent part of DDSCAT.Sa is the single subroutine 
TIMEIT. Several versions of TIMEIT are provided: 

• timeit_sun. f uses the SunOS system call etime 

• timeit_convex. f uses the Convex OS system call etime 

• timeit_cray . f uses the system call second 

• timeitJip.f uses the HP- AUX system calls sysconf and times 

• timeit_ibm6000 . f uses the AIX system call mclock 
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• timeit_osf . f uses the DEC OSF system call etime 

• timeit_sgi . f uses the IRIX system call etime 

• timeit_vms.f uses the VMS system calls LIB$INIT_TIMER and LIB$SHDW_TIMER 

• timeit_titaii.f uses the system call cputim 

• timeit_null . f is a dummy routine which provides no timing information, but can be 
used under any operating system. 

You must compile and link one of the timeit_a;a;a;. f routines, possibly after modifying it to work 
on your system; timeit_null . f is the easiest choice, but it will return no timing information.^ 

6.3 Leaving the netCDF Capability Disabled 

The Makefile supplied with the distribution of DDSCAT.Sa is set up to link to a "dummy" 
subroutine dummywritenet . f instead of subroutine writenet . f , in order to minimize problems 
during initial compilation and linking. The "dummy" routine has no functionality, other than 
bailing out with a fatal error message if the user makes the mistake of trying to specify one of 
the netCDF options (ALLCDF or ORICDF). First-time users of DDSCAT.Sa should not try to 
use the netCDF option - simply skip this section, and specify option NDTCDF in ddscat.par. 



After successfully using DDSCAT.Sa you can return to §3.4 to try to enable the netCDF 
capability. 

6.4 Enabling the netCDF Capability 

Subroutine WRITENET (file writenet. f) provides the capability to output binary data in the 



netCDF standard format (see §10.3). In order to use this routine (instead of dummywritenet . f ), 

it is necessary to take the following steps:^ 

1. Have netCDF already installed on your system (check with your system administrator). 

2. Find out where netcdf . inc is located and edit the include statement in writenet . f to 
specify the correct path to netcdf . inc. 

3. Find out where the libnetcdf .a library is located, and edit the Makefile so that the 
variable LIBNETCDF specifies the correct path to this library. 

4. Edit the Makefile to "comment out" (with a # symbol in column 1) the line writenet = 
dummywritenet and "uncomment" (remove the # symbol) the line writenet = writenet 
so that writenet .f will be compiled instead of the dummy routine dummywritenet .f . 

6.5 Compiling and Linking... 

After suitably editing the Makefile (while still positioned in DDA/src) simply type^ 
make ddscat 

which should create an executable file DDA/src/ddscat . The resulting executable will not have 



^ Non-Unix sites: The VMS-compatible version of TIMEIT is included in the file SRC9 . FOR. For non-VMS sites, you 
will need to replace this version of TIMEIT with one of the other versions, which are to be found in the file MISC. FOR. 
If you are having trouble getting this to work, choose the "dummy" version of TIMEIT from MISC . FDR: this will return 
no timing information, but is platform-independent. 

" Non-UNIX sites: You need to replace the "dummy" version of SUBROUTINE WRITENET in SRC 1. FOR with the 
version provided in MISC. FOR. You will also need to consult your systems administrator to verify that the netCDF 
library has already been i nsta lled on your system, and to find out how to link to the library routines. 

^ Non-Unix sites: see §0.6 
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netCDF capability, and will contain timing instructions compatible with the Solaris 1.x or 2.x 
operating systems on Sun computers, as well as several other versions of Unix. To add netCDF 
capability, see §3.4. To replace the Sun-compatible timing routine with another, see §3.2. 



6.6 Installation on Non-Unix Systems 

DDSCAT.5a is written in standard Fortran-77 plus the DO ... ENDDO extension which ap- 
pears to be supported by all current Fortran-77 compatible compilers. It is possible to run 
DDSCAT on non-Unix systems. If the Unix "make" utility is not available, here in brief is 
what needs to be accomplished: 

All of the necessary Fortran code to compile and link DDSCAT. 5a is included in the 
following files: SRCO.FOR, SRCl.FOR, SRC2.F0R, SRC3.F0R, SRC4.F0R, SRC5.F0R, SRC6.F0R, 
SRC7.F0R, SRC8.F0R, SRC9.F0R, CGCOMMON . FOR, GPFA.FOR, LAPACK.FOR, and PIM.FOR. There 
is an additional file MI SC. FOR, but this is not needed for "basic" use of the code (see below). 

The main program DDSCAT is found in SRCO.FOR. It calls a number of subroutines, which 
are included in the other * . FOR files. 

Three of the subroutines exist in more than one version. The "default" version of each is 
located in the file SRC 1 .FOR: 

• Select the version of SUBROUTINE WRITENET from SRCl.FOR (do not use the version in 
MISC. FOR). The resulting code will not support netCDF capability. (If netCDF capa- 
bility is required, you will need to install netCDF libraries on your system - see §|6^). 

• Select the version of SUBROUTINE C3DFFT from SRCl.FOR (do not use the version in 
MISC. FOR) 

• There is one version of SUBROUTINE TIMEIT included in SRCl.FOR, and a number of ad- 
ditional versions in MISC. FOR. Use the version from SRCl .FOR, which begins 

SUBROUTINE TIMEIT (CMSGTM.DTIME) 

C 

C timeit_null 
C 

C This version of timeit is a dummy which does not provide any 
C timing information. 

This version does not use any system calls, and therefore the code should compile and 
link without problems; however, when you run the code, it will not report any timing 
information reporting how much time was spent on different parts of the calculation. 
If you wish to obtain timing information, you will need to find out what system calls are 
supported by your operating system. You can look at the other versions of SUBROUTINE 
TIMEIT in MISC. FOR to see how this has been done under VMS and various version of 
Unix. 

Once you have selected the appropriate versions of WRITENET, C3DFFT, and TIMEIT, you can sim- 
ply compile and link just as you would with any other Fortran code with a number of modules. 
You should end up with an executable with a (system-dependent) name like DDSCAT.EXE. 

In addition to program DDSCAT, we provide two other Fortran 77 programs which may be 
useful. Program CALLTARGET can be used to call the target generation routines which may be 
helpful for testing purposes. Program TSTFFT is useful for comparing speeds of different FFT 
options. 
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7 Moving the Executable 



Now reposition yourself into the directory DDA (e.g., type cd . . ), and copy the executable from 
src/ddscat to the DDA directory by typing 

cp src/ddscat ddscat 
This should copy the file DDA/src/ddscat to DDA/ddscat . Similarly, copy the sample parameter 
file ddscat .par and the file diel .tab to the DDA directory by typing 

cp doc/ddscat .par ddscat. par 

cp doc/diel.tab diel. tab 



8 The Parameter File ddscat. par 

The directory DDA should now contain a sample file ddscat .par which provides parameters to 
the program ddscat. As provided (see AppendixQ), the file ddscat. par is set up to calculate 
scattering by a 8x6x4 rectangular array of 192 dipoles, with an effective radius ~ l//m, at 
a wavelength of 6.2832/im (for a "size parameter" 27raeff/A = 1). 

The dielectric function of the target material is provided in the file diel. tab, which is a 
sample file in which the refractive index is set to m = 1.33 + O.Oli at all wavelengths; the name 
of this file is provided to ddscat by the parameter file ddscat. par. 

The sample parameter file as supplied calls for the new GPFA FFT routine (GPFAFT) of 
Temperton (1992) to be employed and the PBCGST iterative method to be used for solving the 
system of linear equations. (See section §|l^ and § p^ for discussion of choice of FFT algorithm 
and choice of equation-solving algorithm.) 

The sample parameter file specifies (via option LATTDR) that the "Lattice Dispersion Rela- 
tion" of Draine and Goodman (1993) be employed to determine the dipole polarizabilities. See 
§ pl] for discussion of other options. 

The sample ddscat. par file specifies that the calculations be done for a single wavelength 
(6.2832/j,m) and a single effective radius (1.00/im). Note that in DDSCAT. 5a the "effective 
radius" Oeff is the radius of a sphere of equal volume - i.e., a sphere of volume iVd^ , where d 
is the lattice spacing and N is the number of occupied (i.e., non- vacuum) lattice sites in the 
target. Thus the effective radius a^s = (37V/47r)^/^d . 

The incident radiation is always assumed to propagate along the x axis in the "Lab Frame" . 
The sample ddscat . par file specifies incident polarization state egi to be along the y axis (and 
consequently polarization state eo2 will automatically be taken to be along the z axis). I0RTH=2 
in ddscat .par calls for calculations to be carried out for both incident polarization states (egi 
and eo2 - see ^ ). 

The target is assumed to have two vectors ai and 3.2 embedded in it; 3.2 is perpendicular to 
3.1. In the case of the 8x6x4 rectangular array of the sample calculation, the vector ai is along 
the "long" axis of the target, and the vector a2 is along the "intermediate" axis. The target 
orientation in the Lab Frame is set by three angles: /3, 6, and $, defined and discussed below 
in § [r7| . Briefly, the polar angles and $ specify the direction of ai in the Lab Frame. The 
target is assumed to be rotated around ai by an angle (3. The sample ddscat .par file specifies 
(3 — Q and $ = (see lines in ddscat. par specifying variables BETA and PHI), and calls for 
three values of the angle Q (see line in ddscat. par specifying variable THETA). DDSCAT. 5a 
chooses 8 values uniformly spaced in cos0; thus, asking for three values of 8 between and 
90° yields 8 = 0, 60°, and 90°. 

Appendix ^ provides a detailed description of the file ddscat. par 

Note that the structure of ddscat .par depends on the version of DDSCAT being used. Make sure you update 
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9 Running DDSCAT.5a Using the Sample ddscat .par File 



To execute the program on a UNIX system (running cither sh or csh), simply type 

ddscat >& ddscat. out & 
which will redirect the "standard output" to the file ddscat. out, and run the calculation in 
the background. The sample calculation (8x6x4=192 dipole target, 3 orientations, two incident 
polarizations, with scattering calculated for 14 distinct scattering directions), requires 1.0 cpu 
seconds on a Sun Ultra- 170 workstation. 

10 Output Files 
10.1 ASCII files 

If you run DDSCAT using the command 
ddscat >& ddscat. out & 

you will have various types of ASCII files when the computation is complete: 

• a file ddscat . out; 

• a file mtable; 

• a file qtable; 

• a file qtable2; 

• files vxxryyorl . avg (one, wOOrOOori . avg, for the sample calculation); 

• if ddscat . par specified IWRKSC=1, there will also be files vixxryykzzz. sea (3 for the sample 
calculation: wOOrOOkOOO . sea, wOOrOOkOOl . sea, w00r00k002 . sea). 

The file ddscat. out will contain any error messages generated as well as a running report on 
the progress of the calculation, including creation of the target dipole array. During the iterative 
calculations, Qoxt, Qabs, and Qpha are printed after each iteration; you will be able to judge the 
degree to which convergence has been achieved. Unless TIMEIT has been disabled, there will 
also be timing information. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Qext, Qabs, 
Qsca, g(l) = (cos(0s)), and Qbk- Here Qcxt, Qabs, and Qsca are the extinction, absorption, and 
scattering cross sections divided by Tra^g. Qbk is the differential cross section for backscattering 
(area per sr) divided by Tra^g. 

The file qtable2 contains a summary of the orientationally-averaged values of Qpha, Qpoi, 
and Qcpoi- Here Qpha is the "phase shift" cross section divided by na^g (see definition in 
Draine 1988). Qpoi is the "polarization efficiency factor", equal to the difference between Qext 
for the two orthogonal polarization states. We define a "circular polarization efficiency factor" 
Qcpoi = QpoiQpha, since an optically-thin medium with a small twist in the alignment direction 
will produce circular polarization in initially unpolarized light in proportion to Qcpoi- 

For each wavelength and size, DDSCAT. 5a produces a file with a name of the form 
wa;2xi/i/ori . avg, where index xx (=00, 01, 02....) designates the wavelength and index yy (=00, 
01, 02...) designates the "radius" ; this file contains Q values and scattering information averaged 
over however many target orientations have been specified (see §17|. The file wOOrOOori . avg 
produced by the sample calculation is provided below in Appendix 

old parameter files before using them with DDSCAT. 5a ! 
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In addition, if ddscat . par has specified IWRKSC=1 (as for the sample calculation), DDSCAT.5a 
will generate files with names of the form vixxryifs.zzz.BMg, where xx and yy are as before, and 
index zzz =(000,001,002...) designates the target orientation; these files contain Q values and 
scattering information for each of the target orientations. The structure of each of these files 
is very similar to that of the vxxryyorl . avg files. Because these files may not be of particular 
interest, and take up disk space, you may choose to set IWRKSC=0 in future work. However, it 
is suggested that you run the sample calculation with IWRKSC=1. 

The sample ddscat. par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 tar- 
get size, and averaging over 3 target orientations. Running DDSCAT. 5a with the sample 
ddscat .par file will therefore generate files wOOrOOkOOO . sea, wOOrOOkOOl . sea, and w00r00k002 . sea 
. To understand the information contained in one of these files, please consult Appendix 
which contains an example of the file wOOrOOkOOO . sea produced in the sample calculation. 



10.2 Binary Option 

It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete infor- 
mation, including header and data sections. This is accomplished by specifying either ALLBIN 
or ORIBIN in ddseat .par . 

Subroutine writebin.f provides an example of how this can be done. The "header" sec- 
tion contains dimensioning and other variables which do not change with wavelength, particle 
geometry, and target orientation. The header section contains data defining the particle shape, 
wavelengths, particle sizes, and target orientations. If ALLBIN has been specified, the "data" 
section contains, for each orientation, Mueller matrix results for each scattering direction. The 
data output is limited to actual dimensions of arrays; e.g. nseat,4,4 elements of Mueller 
matrix are written rather than mxseat,4,4. This is an important consideration when writing 
postprocessing codes. 

A skeletal example of a postprocessing code was written by us and is provided in sub- 
directory DDA/IDL. If you do plan to use the Interactive Data Language (IDL) for postpro- 
cessing, you may consider using the netCDF binary file option which offers substantial ad- 
vantages over the FORTRAN unformatted write. More information about IDL is available at 
http : //www . rsine . com/ idl. Unfortunately IDL requires a license and hence is not distributed 
with DDSCAT. 



10.3 Machine-Independent Binary File Option: netCDF 

The "unformatted" binary file is compact, fairly complete, and (in comparison to ASCII output 
files) is efficiently written from FORTRAN. However, binary files are not compatible between 
different computer architectures if written by "unformatted" WRITE from FORTRAN. Thus, you 
have to process the data on the same computer architecture that was used for the DDSCAT 
calculations. We provide an option of using netCDF with DDSCAT. The netCDF library defines 
a machine-independent format for representing scientific data. Together, the interface, library, 
and format support the creation, access, and sharing of scientific data. For more information 



see http : //www . unidata. uear . edu/paekages/netedf 



Several major graphics packages (for example IDL) have adopted netCDF as a standard 
for data transport. For this reason, and because of strong and free support of netCDF over 
the network by UNIDATA, we have implemented a netCDF interface in DDSCAT. This may 
become the method of choice for binary file storage of output from DDSCAT. 

After the initial "learning curve" netCDF offers many advantages: 

• It is easy to examine the contents of the file (using tools such as nedump) . 
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• Binary files are portable - they can be created on a supercomputer and processed on a 
workstation. 

• Major graphics packages now offer netCDF interfaces. 

• Data input and output are an order of magnitude faster than for ASCII files. 

• Output data files are compact. 
The disadvantages include: 

• Need to have netCDF installed on your system. 

• Lack of support of complex numbers. 

• Nontrivial learning curve for using netCDF. 

• Lack of portability of netCDF libraries. 

The public-domain netCDF software is available for many operating systems from 



tittp: //www. unidata. ucar.edu/packages/netcdf . The steps necessary for enabling the 



netCDF capability in DDSCAT.Sa are listed above in ^6.4. Once these have been success- 
fully accomplished, and you are ready to run DDSCAT to produce netCDF output, you must 
choose either the ALLCDF or ORICDF option in ddscat .par; ALLCDF will result in a file which will 
include the Mueller matrix for every wavelength, particle size, and orientation, whereas DRICDF 
will result in a file limited to just the orientational averages for each wavelength and target size. 



11 Dipole Polarizabilities 

Option LATTDR specifies that the "Lattice Dispersion Relation" of Draine and Goodman (1993) 
be employed to determine the dipole polarizabilities.; other possible choices are DRAI88 (pre- 
scription used by Draine 1988) and GDBR88 (prescription used by Goedecke & O'Brien 1988 
and Hage & Greenberg 1990). In the limit |m|fc(i <C 1 (where k = 27r/A is the wave vector 
and d is the lattice spacing) all three options converge to the same limit; for |TO|/c(i > 0.1 there 
are significant differences among them. Consult the paper by Draine & Goodman (1993) for 
discussion and comparison of these three prescriptions. Option LATTDR is recommended for 
general use, based upon the tests presented by Draine & Goodman (1993). 



12 Dielectric Functions 

In order to assign the appropriate dipole polarizabilities, DDSCAT.Sa must be given the 
dielectric constant of the material (or materials) of which the target of interest is composed. 
Unless the user wishes to use the dielectric function of either solid or liquid H2O (see below), his 
information is supplied to DDSCAT through a table (or tables), read by subroutine DIELEC 
in file dielec.f , and providing either the complex refractive index m = n + ik or complex 
dielectric function e = ei + 162 as a function of wavelength A. Since m = e^/^, or e = m^, the 
user must supply either to or e, but not both. 

The table formatting is intended to be quite flexible. The first line of the table consists of 
text, up to 80 characters of which will be read and included in the output to identify the choice 
of dielectric function. (For the sample problem, it consists of simply the statement m = 1 . 33 
+ O.Oli.) The second line consists of 5 integers; either the second and third or the fourth and 
fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(TO) is stored in. 
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• The third integer specifies which column Im(m) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 



• The fifth integer specifies which column Im(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and Im(e) from the file; 
if the fourth and fifth integers are zeros, then Re(m) and Im(m) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There 
must be at least 3 lines of data: even if e is required at only one wavelength, please supply two 
additional "dummy" wavelength entries in the table so that the interpolation apparatus will 
not be confused. 

In the event that the user wishes to compute scattering by targets with the refractive index of 
either solid or liquid H2O, we have included two "built-in" dielectric functions. If H2DICE is spec- 
ified on line 10 of ddscat .par, DDSCAT will use the dielectric function of H2O ice at T=250K 
as compiled by Warren (1984). If H20Liq is specified on fine 10 of ddscat .par, DDSCAT will 
use the dielectric function for liquid II2O at T=280K using subroutine REFWAT written by Eric 
A. Smith. For more information sec tLttp://atol.ucsd.edu/^pflatau/scatlib/refr .html . 



13 Choice of FFT Algorithm 

One major change in going from DDSCAT. 4b to 4c and 5a was modification of the code 
to permit use of the GPFA FFT algorithm developed by Dr. Clive Temperton. DDSCAT 
continues to offer both the Brenner code as well as the "old" Temperton code as alternative 
FFT implementations. The "old" Temperton code requires about 11% more memory than 
either the Brenner or GPFA codes (to use the "old" Temperton algorithm, DDSCAT . f must be 
compiled with MXMEM=:1 rather than MXMEM^O). 

To help persuade the user that the GPFA code is an important step forward, we provide a 
program TSTFFT to compare the performance of different 3-D FFT implementations. To compile, 
link, and run this program on a Unix system,^ position yourself in the DDA/ src directory and 
type 

make tstfft 
tstf ft 

Output results will be written into a file res.dat. Here is a copy of the res.dat file created 
when run on a Sun Ultrasparc 170: 

CPU time (sec) for different 3-D FFT methods 
Machine= Sun Ultrasparc-170 
parameter LVR = 64 

LVR="length of vector registers" in gpf a2f ,gpf a3f ,gpf a5f 









Brenner 


Temperton 


Temperton 


NX 


NY 


NZ 


(FOURX) 


(Old) 


(GPFA) 


2 


2 


2 


0.000068 


0.000065 


0.000262 


3 


3 


3 


0.000200 


0.000077 


0.000096 


4 


4 


4 


0.000064 


0.000086 


0.000111 


5 


5 


5 


0.000651 


0.000292 


0.000137 


6 


6 


6 


0.001171 


0.000170 


0.000223 


8 


8 


8 


0.000350 


0.000250 


0.000323 



Non-Unix systems; the TSTFFT Fortran source code is in the file MISC. FOR. 
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9 


9 


9 


0, 


.004654 


0, 


.000300 


0, 


.000505 


10 


10 


10 


0, 


.005218 


0, 


.000457 


0, 


.000575 


12 


12 


12 


0, 


.008664 


0, 


.000814 


0, 


.000751 


15 


15 


15 


0, 


.023387 


0, 


.001583 


0, 


.001798 


16 


16 


16 


0, 


.003137 


0, 


.001359 


0, 


.002386 


18 


18 


18 


0, 


.042196 


0, 


.003908 


0, 


.003479 


20 


20 


20 


0, 


.043754 


0, 


.003881 


0, 


.004010 


24 


24 


24 


0, 


.075672 


0, 


.010722 


0, 


.007198 


25 


25 


25 


0, 


.103492 


0, 


.008142 


0, 


.011232 


27 


27 


27 


0, 


.160173 


0, 


.011126 


0, 


.012293 


30 


30 


30 


0, 


.185177 


0, 


.017111 


0, 


.014631 


32 


32 


32 


0, 


.042211 


0, 


.034032 


0, 


.023245 


36 


36 


36 


0, 


.322581 


0, 


.071013 


0, 


.027502 


40 


40 


40 


0, 


.379551 


0, 


.133391 


0, 


.044182 


45 


45 


45 


0, 


.797107 


0, 


.199537 


0, 


.069740 


48 


48 


48 


0, 


.668408 


0, 


.255849 


0, 


.094520 


50 


50 


50 


0, 


.983722 


0, 


.287208 


0, 


.115833 


54 


54 


54 


1, 


.530319 


0, 


.427839 


0, 


. 142388 


60 


60 


60 


1, 


.804154 


0, 


.607189 


0, 


.181221 


64 


64 


64 


1, 


.005013 


0, 


.806311 


0, 


.342762 


72 


72 


72 


3, 


.812469 


1, 


.219364 


0, 


.362328 


75 


75 


75 


5, 


.364067 


1, 


.234528 


0, 


.437496 


80 


80 


80 


4, 


.509048 


1, 


.576118 


0, 


.563839 


81 


81 


81 


8, 


.017456 


1, 


.798195 


0, 


.560326 


90 


90 


90 


10, 


.201881 


2, 


.455399 


0, 


.779869 



It is clear that the GPFA code is generally much faster than the Brenner FFT (by a factor 
of 13 for the 90x90x90 case) and significantly faster than the "old" Temperton FFT code (by 
a factor of 3 for the 90x90x90 case), although for some cases the differences are not large (e.g., 
27x27x27). Since the GPFA code is memory-efhcient as well, it appears to be the method of 
choice on scalar machines. It appears to also be best on vector machines. 

The GPFA code contains a parameter LVR which is set in data statements in the routines 
gpf a2f , gpfa3f , and gpfaSf . LVR is supposed to be optimized to correspond to the "length 
of a vector register" on vector machines. As delivered, this parameter is set to 64, which is 
supposed to be appropriate for Grays other than the G90 (for the G90, 128 is supposed to be 
preferable) .0 The value of LVR is not critical for scalar machines, as long as it is fairly large. 
We found little difference between LVR=64 and 128 on a Sparc 10/51 and on an Ultrasparc 
170. You may wish to experiment with different LVR values on your computer architecture. To 
change LVR, you need to edit gpfa.f and change the three data statements where LVR is set. 

The choice of FFT implementation is obtained by specifying one of: 

• GPFAFT to use the GPFA algorithm (Temperton 1992) - this is recommended!; 

• TMPRTN to obtain the "old" Temperton (1983) implementation; 

• BRENNR to obtain the FOURX implementation of the Brenner (1969) algorithm; 

• CONVEX to use the native FFT routine on a Convex system 

* In place of "preferable" users are encouraged to read "necessary" - there is some basis for fearing that results 
computed on a C90 with LVR other than 128 run the risk of being incorrect! Please use LVR=128 if running on a C90! 
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14 Choice of Iterative Algorithm 



As discussed elsewhere (e.g., Drainc 1988), the problem of electromagnetic scattering of an 
incident wave Ej„c by an array of N point dipoles can be cast in the form 

AP = E (17) 

where E is a 3iV-dimensional (complex) vector of the incident electric field 'Emc at the N lattice 
sites, P is a 3iV-dimensional (complex) vector of the (unknown) dipole polarizations, and A is 
a 3N X 3iV complex matrix. 

Because 3N is a large number, direct methods for solving this system of equations for 
the unknown vector P are impractical, but iterative methods are useful: we begin with a 
guess (typically, P = 0) for the unknown polarization vector, and then iteratively improve the 
estimate for P until equation (|l^) is solved to some error criterion. The error tolerance may be 
specified as 

|AtAP-AtE| 

' 

where A^ is the Hermitian conjugate of A [(A^)ij = (Aji)*], and h is the error tolerance. We 
typically use h — 10^^ in order to satisfy eq.(p^ to high accuracy. The error tolerance h can 
be specified by the user (see Appendix ^) . 

A major change in going from DDSCAT.4b to 5a is the implementation of several different 
algorithms for iterative solution of the system of complex linear equations. DDSCAT.Sa 
has been modified to permit solution algorithms to be treated in a fairly "modular" fashion, 
facilitating the testing of different algorithms. Many algorithms were compared by Flatau 
(1997)^; two of them performed well and are made available to the user in DDSCAT.Sa. The 
choice of algorithm is made by specifying one of the options: 

• PBCGST ~ Preconditioned BiConjugate Gradient with STabilitization method from the 
Parallel Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), 
as coded in the Complex Conjugate Gradient package (CCGPACK) created by P.J. 
Flatau. This is the algorithm discussed by Draine (1988) and used in previous versions of 
DDSCAT. 

Both methods work well. Our experience suggests that PBCGST is often faster than PETRKP, by 
perhaps a factor of two. We therefore recommend it, but include the PETRKP method as an 
alternative. 

The Parallel Iterative Methods (PIM) by Rudnci Dias da Cunha (rddSukc . ac . uk) and Tim 
Hopkins (trh@ukc.ac.uk) is a collection of Fortran 77 routines designed to solve systems of 
linear equations on parallel and scalar computers using a variety of iterative methods (available 
at 



tittp : //www. mat .uf rgs .br/pim-e .html). PIM offers a number of iterative methods, including 

• Conjugate-Gradients, CG (Hestenes 1952), 

• Bi-Conjugate-Gradients, BIGG (Fletcher 1976), 

• Conjugate-Gradients squared, CCS (Sonneveld 1989), 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991), 

• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992) 



A postscript copy of this report - file cg.ps - is distributed with the DDSCAT.Sa documentation. 
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• the restarted generalized minimal residual, RGMRES (Saad 1986), 

• the restarted generalized conjugate residual, RGCR (Eisenstat 1983), 

• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955), 

• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996), 

• transpose-free quasi-minimal residual, TFQMR (Freund 1992), 

• the Chebyshev acceleration, CHEBYSHEV (Young 1981). 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can 
be called directly via ddscat . par. It is possible (and was done) to add other options by changing 
the code in getf ml . f . A helpful introduction to conjugate gradient methods is provided by the 
report "Conjugate Gradient Method Without Agonizing Pain" by Jonathan R. Shewchuk, avail- 



able as a postscript file: ftp : //REPORTS . ADM . CS . CMU . EDU/usr0/anon/1994/CMU-CS-94-125 . ps 



15 Calculation of Radiative Force and Torque 

In addition to solving the scattering problem for a dipole array, DDSCAT. 5a can compute 
the three-dimensional force Frad and torque Fjad exerted on this array by the incident and 
scattered radiation fields. This calculation is carried out, after solving the scattering problem, 
provided DDTDRQ has been specified in ddscat. par. For each incident polarization mode, the 
results are given in terms of dimensionless efficiency vectors Qpr and Qr, defined by 

Q^^JEl^ , (20) 

Tra^gUrad 

where Frad and Fiad are the time-averaged force and torque on the dipole array, k = 2tt/X 
is the wavenumber in vacuo, and Uiad — Eq/Stt is the time-averaged energy density for an 
incident plane wave with amplitude Eq cos{ujt + (p). The radiation pressure efficiency vector can 
be written 

Qpr — Qcxtk QscaS : (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean 
direction of propagation of the scattered radiation: 

1 /•^^ dCsca(n,k) . 

where d£l is the element of solid angle in scattering direction n, and dCsca./dfl is the differential 
scattering cross section. 

Equations for the evaluation of the radiative force and torque are derived by Draine & Wein- 
gartner (1996). It is important to note that evaluation of Qpr and Qr involves averaging over 
scattering directions to evaluate the linear and angular momentum transport by the scattered 
wave. This evaluation requires appropriate choices of the parameters ICTHM and IPHM - see g23. 
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16 Memory Requirements 



Since Fortran-77 does not allow for dynamic memory allocation, the executable has compiled 
into it the dimensions for a number of arrays; these array dimensions limit the size of the dipole 
array which the code can handle. The source code supplied to you (which can be used to run the 
sample calculation with 192 dipoles) is restricted to problems with targets which are maximum 
extent of 8 lattice spacings along the x-, y-, and z-directions (MXNX=8 ,MXNY=8 ,MXNZ=8; i.e, the 
target must fit within an 8x8x8=512 cube) and involve at most 9 different dielectric functions 
(MXCDMP=9). With this dimensioning, the executable requires about 1.3 MB of memory to run 
on an Ultrasparc system; memory requirements on other hardware and with other compilers 
should be similar. To run larger problems, you will need to edit the code to change PARAMETER 
values and recompile. 

All of the dimensioning takes place in the main program DDSCAT - this should be the only rou- 
tine which it is necessary to recompile. All of the dimensional parameters are set in PARAMETER 
statements appearing before the array declarations. You need simply edit the parameter state- 
ments. Remember, of course, that the amount of memory allocated by the code when it runs 
will depend upon these dimensioning parameters, so do not set them to ridiculously large val- 
ues! The parameters MXNX, MXNY, MXNZ specify the maximum extent of the target in the x-, y-, 
or z-directions. Set the parameter MXMEM=0 or 1 depending on whether workspace required by 
Temperton's FFT algorithm is to be reserved. The parameter MXCOMP specifies the maximum 
number of different dielectric functions which the code can handle at one time. The comment 
statements in the code supply all the information you should need to change these parameters. 

The memory requirement for DDSCAT. 5a (with the netCDF option disabled) is approx- 
imately (1059 + 0.623MXNX X MXNY X MXNZ) kbytes when compiled with MXMEM=0, or (1059 -f 
0.686MXNXX MXNY X MXNZ) kbytes when compiled with MXMEM=1. Thus with MXMEM=0 a 32x32x32 
calculation requires 21.0 MBytes, while a 48x48x48 calculation requires 68.4 MBytes. These 
values are for an Ultrasparc system using the Sun FORTRAN compiler. 

17 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x 
direction. In ddscat .par one specifies the first polarization state egi (which obviously must lie 
in the y,z plane in the LF); DDSCAT automatically constructs a second polarization state 
eo2 = X X Gq]^ orthogonal to epi (here x is the unit vector in the +x direction of the LF. For 
purposes of discussion we will always let unit vectors x, y, z = x x y be the three coordinate 
axes of the LF. Users will often find it convenient to let polarization vectors Bqi = y, eo2 = z 
(although this is not mandatory - see §|2l|). 

Recall that definition of a target involves specifying two unit vectors, ai and 3.2, which are 
imagined to be "frozen" into the target. We require a2 to be orthogonal to ai. Therefore we 
may define a "Target Frame" (TF) defined by the three unit vectors ai , a2 , and as = ai x a2 . 

For example, when DDSCAT creates a 8x6x4 rectangular solid, it fixes ai to be along the 
longest dimension of the solid, and 3.2 to be along the next-longest dimension. 

Orientation of the target relative to the incident radiation can in principle be determined 
two ways: 

1. specifying the direction of ai and a2 in the LF, or 

2. specifying the directions of x (incidence direction) and y in the TF. 

DDSCAT. 5a uses method 1.: the angles 0, $, and /3 are specified in the file ddscat .par. 
The target is oriented such that the polar angles 9 and $ specify the direction of ai relative to 
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Figure 5: Target orientation in the Lab Frame, x is the direction of propagation of the incident 
radiation, and y is the direction of the "real" component of the first incident polarization mode. 
In this coordinate system, the orientation of target axis ai is specified by angles and With 
target axis ai fixed, the orientation of target axis a2 is then determined by angle (3 specifying 
rotation of the target around ai. When /3 = 0, a2 lies in the ai,x plane. 
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the incident direction x, where the x,y plane has $ = 0. Once the direction of ai is specified, 
the angle /? than specifies how the target is to rotated around the axis ai to fully specify its 
orientation. A more extended and precise explanation follows: 



17.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 9, $, and (3, to specify the directions of unit vectors ai and a2 in 
the LF (see Fig. |). 

Q is the angle between ai and x. 

When <I> = 0, ai will lie in the x, y plane. When $ is nonzero, it will refer to the rotation 
of ai around x: e.g., $ = 90° puts ai in the x, z plane. 

When P — 0, a2 will lie in the x, ai plane, in such a way that when 8 = and $ = 0, 3.2 is 
in the y direction: e.g, 8 = 90°, $ = 0, /3 = has ai = y and a2 = — x. Nonzero f3 introduces 
an additional rotation of a2 around ai: e.g., 8 = 90°, $ = 0, /3 = 90° has ai = y and a2 = z. 

Mathematically : 

ai = X cos 8 + y sin 8 cos $ + z sin 8 sin $ (23) 
a2 ~ — xsin8cos/3 + y[cos8cos/3cos$ — sin/3sin<l>] 

+z [cos 8 cos /? sin $ + sin /? cos $] (24) 
as — xsin8sin/3 — y[cos8sin/3cos<i> + cos/Jsin^] 

— z[cos8sin/3sin<i> — cos/3cos$] (25) 



or, equivalently: 



X = ai cos 8 — a2 sin 8 cos /3 + as sin 8 sin /3 (26) 
y = ai sin 8 cos $ + a2 [cos 8 cos /3 cos $ — sin (3 sin $] 

-a3[cos8sin/3cos$ + cos/3sin<I>] (27) 
z = ai sin 8 sin $ + a2 [cos 8 cos (3 sin $ + sin (3 cos $] 

—as [cos 8 sin /3 sin $ — cos (3 cos $] (28) 



17.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that x (the 
direction of propagation of the radiation) and y (usually the first polarization direction) and z 
(= X X y) refer to certain directions in the TF. Given the definitions of the LF and TF above, 
this is simply an exercise in coordinate transformation. For example, one might wish to have 
the incident radiation propagating along the (1,1,1) direction in the TF (example 14 below). 
Here we provide some selected examples: 

1. X = ai, y = a2, z = as : 8 = 0, $ + /? = 

2. X = ai, y = as, z = -a2 : 8 = 0,$ + /? = 90° 

3. X = a2, y = ai, z = -as : 8 = 90°, f3 = 180°, $ = 

4. X = a2, y = as, z = ai : 8 = 90°, (3 = 180°, $ = 90° 

5. X = as, y = ai, z = a2 : 8 = 90°, (3 = -90°, $ = 

6. X = as, y = aa, z = -ai : 8 = 90°, f3 = -90°, $ = -90° 

7. X = -ai, y = a2, z = -as : 8 = 180°, /3 + $ = 180° 

8. X = -ai, y = as, z = £12 : 8 = 180°, /? + $ = 90° 
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9. 


X = 


-a2, y = ai, z as : e = 


90°, /3 = 0, $ = 




10. 


X = 


-a2, y = as, z = -ai : 9 


= 90°, /3 = 0, $ = - 


-90° 


11. 


X = 


-as, y = ai, z = -as : 9 


= 90°, /3 = -90°, $ 


= 


12. 


X = 


-as, y = as, z = ai : 9 = 


90°, /3 = -90°, $ = 


: 90° 


13. 


ic — 


(ai +as)/\/2, y = as, z = 


(ai - as)/V2: 6 = 


45°, /3 = 180°, $ 90° 


14. 


X = 

/3 = 


(ai + as + a3)/y3, y = 
135°, $ = 30°. 


(ai - a2)/V2, z = 


(ai + as - 2as)/y6: 6 



17.3 Sampling in 0, $, and (i 

The present version, DDSCAT.Sa, chooses the angles /3, 9, and $ to sample the inter- 
vals (BET AMI, BET AMX), (THETMI ,THETMX) , (PHIMIN.PHIMAX), where BETAMI, BETAMX, THETMI, 
THETMX, PHIMIN, PHIMAX are specified in ddscat .par . The prescription for choosing the angles 
is to: 

• uniformly sample in /3; 

• uniformly sample in <&; 

• uniformly sample in cos 9 . 

This prescription is appropriate for random orientation of the target, within the specified limits 
of /3, and 9. 

Note that when DDSCAT.Sa chooses angles it handles /3 and $ differently from 9.0 The 
range for /? is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if 
you take BETAMI=0, BETAMX=90, NBETA=2 you will get /3 = 22.5° and 67.5°. Similarly, if you 
take PHIMIN=0, PHIMAX=180, NPHI=2 you wiU get $ = 45° and 135°. 

Sampling in 9 is done quite differently from sampling in /3 and <i>. First, as already mentioned 
above, DDSCAT.Sa samples uniformly in cos 9, not 9. Secondly, the sampling depends on 
whether NTHETA is even or odd. 

• If NTHETA is odd, then the values of 9 selected include the extreme values THETMI and 
THETMX; thus, THETMI=0, THETMX=90, NTHETA=3 will give you 9 = 0,60°, 90°. 

• If NTHETA is even, then the range of cos 9 will be divided into NTHETA intervals, and the 
midpoint of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give 
you 9 41.41° and 75.52° [cos9 = 0.25 and 0.75]. 

The reason for this is that if odd NTHETA is specified, then the "integration" over cos 9 is 
performed using Simpson's rule for greater accuracy. If even NTHETA is specified, then the 
integration over cos 9 is performed by simply taking the average of the results for the different 
9 values. 

If averaging over orientations is desired, it is recommended that the user specify an odd 
value of NTHETA so that Simpson's rule will be employed. 

18 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How 
to go about this depends on the distribution of orientations which is applicable. 

This is a change from DDSCAT.4a!!. 
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18.1 Randomly-Oriented Targets 



For randomly-oriented targets, we wish to compute the orientational average of a quantity 
Q(/3,e,<i>): 

(Q)^7rT dp dcose / d$Q(/3,e,$) . (29) 

Stt^ Jo J_i Jo 

To compute such averages, all you need to do is edit the file ddscat.par so that DDSCAT 
knows what ranges of the angles /3, Q, and $ arc of interest. For a randomly-oriented target 
with no symmetry, you would need to let /3 run from to 360°, 6 from to 180°, and <i> from 
to 360°. 

For targets with symmetry, on the other hand, the ranges of /3, 0, and $ may be reduced. 
First of all, remember that averaging over $ is relatively "inexpensive" , so when in doubt 
average over to 360°; most of the computational "cost" is associated with the number of 
different values of {(3,Q) which are used. Consider a cube, for example, with axis ai normal to 
one of the cube faces; for this cube /3 need run only from to 90°, since the cube has fourfold 
symmetry for rotations around the axis ai. Furthermore, the angle Q need run only from to 
90°, since the orientation (/3,e,$) is indistinguishable from (/?, 180° - 6, 360° - $). 

For targets with symmetry, the user is encouraged to test the significance of /3,8,$ on targets 
with small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

It is important to remember that DDSCAT. 4b handled even and odd values of NTHETA 
differently - see ^ above! For averaging over random orientations odd values of NTHETA are to 
be preferred, as this will allow use of Simpson's rule in evaluating the "integral" over cos0. 

18.2 Nonrandomly-Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around 
the X axis = direction of propagation of the incident radiation), one may be able to use 
DDSCAT. 5a with appropriate choices of input parameters. More generally, however, you will 
need to modify subroutine ORIENT to generate a list of NBETA values of /?, NTHETA values of 9, 
and NPHI values of $, plus two weighting arrays WGTA(1-NTHETA, 1-NPHI) and WGTB(l-NBETA) . 
Here WGTA gives the weights which should be attached to each (9,$) orientation, and WGTB gives 
the weight to be attached to each P orientation. Thus each orientation of the target is to be 
weighted by the factor WGTAxWGTB. For the case of random orientations, DDSCAT. 5a chooses 
9 values which are uniformly spaced in cos 9, and /3 and $ values which are uniformly spaced, 
and therefore uses uniform weights 

WGTB=1. /NBETA 
When NTHETA is even, DDSCAT sets 

WGTA= 1 . / (NTHETA X NPHI ) 
but when NTHETA is odd, DDSCAT uses Simpson's rule when integrating over 9 and 

WGTA= (1/3 or 4/3 or 2/3)/(NTHETAxNPHl) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly 
oriented cases. If, for example, the orientation is such that for a given $ value only one 9 value 
is possible (this situation might describe ice needles oriented with the long axis perpendicular 
to the vertical in the Earth's atmosphere, illuminated by the Sun at other than the zenith) then 
it is foolish to consider all the combinations of 9 and which the present version of DDSCAT 
is set up to do. We hope to improve this in a future version of DDSCAT. 
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19 Target Generation 



DDSCAT contains routines to generate dipole arrays representing targets of various geometries, 
including spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two 
touching ellipsoids, and three touching ellipsoids. The target type is determined by the variable 
CSHAPE (= ELLIPS, CYLNDR, etc., ...) on line 7 of ddscat .par, and up to 6 parameters (SHPARl, 
SHPAR2, SHPAR3, ...) on line 8. The target geometry is most conveniently described in a 
coordinate system attached to the target which we refer to as the "Target Frame" (TF), so 
in this section only we will let x,y,z be coordinates in the Target Frame. Once the target is 
generated, the orientation of the target in the Lab Frame is accomplished as described in § 1171 . 
Target geometries currently supported include: 

• ELLIPS - Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARl, SHPAR2, SHPAR3 in the x, y, z directions in the TF. (x/SHPARl)^ + 
(y/SHPAR2)2 + (z/SHPARS)^ = d'^/A, where d is the interdipole spacing. 
The target axes are set to ai = (1, 0, 0) and 3.2 — (0, 1, 0) in the TF. 
User must set NCDMP=1 on line 9 of ddscat .par. 

A homogeneous, isotropic sphere is obtained by setting SHPARl = SHPAR2 = SHPAR3 

= diametcr/d. 

• ANIELL Homogeneous, anisotropic ellipsoid. 

SHPARl, SHPAR2, SHPAR3 have same meaning as for ELLIPS. Target axes ai = (1, 0, 0) and 
£12 = (0,1,0) in the TF. It is assumed that the dielectric tensor is diagonal in the TF. 
User must set NCDMP=3 and provide xx, yy, zz elements of the dielectric tensor. 

• CDNELL - Two concentric ellipsoids. 

SHPARl, SHPAR2, SHPAR3 specify the lengths along x, y, z axes (in the TF) of the outer 
ellipsoid; SHPAR4, SHPAR5, SHPAR6 are the lengths, along the x, y, z axes (in the TF) of 
the inner ellipsoid. The "core" within the inner ellipsoid is composed of isotropic material 
1; the "mantle" between inner and outer ellipsoids is composed of isotropic material 2. 
Target axes ai = (1,0,0), a2 = (0,1,0) in TF. User must set NCDMP=2 and provide 
dielectric functions for "core" and "mantle" materials. 

• CYLNDR Homogeneous, isotropic cylinder. 

Length/(i=SHPARl, diameter/d=SHPAR2, with cylinder axis = ai = (1,0,0) and a2 = 
(0, 1, 0) in the TF. User must set NC0MP=1. 

• UNICYL - Homogeneous cylinder with unixial anisotropic dielectric tensor. 

SHPARl, SHPAR2 have same meaning as for CYLNDR. Cyhnder axis = ai = (1,0,0), a2 = 
(0, 1,0). It is assumed that the dielectric tensor e is diagonal in the TF, with tyy = t^z- 
User must set NC0MP=2. Dielectric function 1 is for E || ai (cylinder axis), dielectric 
function 2 is for E _L ai. 

• RCTNGL Homogeneous, isotropic, rectangular solid. 

X, y, z lengths/d = SHPARl, SHPAR2, SHPAR3. Target axes ai = (1,0,0) and a2 = (0, 1,0) 
in the TF. User must set NC0MP=1. 

• HEXGON - Homogeneous, isotropic hexagonal prism. 

SHPARl=Length/(i, SHPAR2 = 2x (hexagon side)/(i = (distance between opposite vertices 
of hexagon)/o? — "hexagon diameter" /d. Target axis ai = (1, 0, 0) in the TF is along the 
prism axis, and target axis 3.2 — (0,1,0) in the TF is normal to one of the rectangular 
faces of the hexagonal prism. User must set NCDMP=1. 

• TETRAH - Homogeneous, isotropic tetrshedron. 

SHPARl=length/c? of one edge. Orientation: one face parallel to y,z plane (in the TF), 
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opposite "vertex" is in +x direction, and one edge is parallel to z axis (in the TF). Target 
axes ai = (1,0,0) [emerging from one vertex] and a2 = (0, 1,0) [emerging from an edge] 
in the TF. User must set NC0MP=1. 

TWOSPH - Two touching homogeneous, isotropic spheroids, with distinct com- 
positions. 

First spheroid has length SHPARl along symmetry axis, diameter SHPAR2 perpendicular 
to symmetry axis. Second spheroid has length SHPAR3 along symmetry axis, diameter 
SHPAR4 perpendicular to symmetry axis. Contact point is on line connecting ccntroids. 
Line connecting centroids is in x direction. Symmetry axis of first spheroid is in y di- 
rection. Symmetry axis of second spheroid is in direction ycos(SHPAR5) + zsin(SHPAR5), 
where y and z arc basis vectors in TF, and SHPAR5 is in degrees. If SHPAR6=0., then target 
axes ai = (1,0,0), a2 = (0, 1,0). If SHPAR6=1., then axes ai and a2 are set to principal 
axes with largest and 2nd largest moments of inertia assuming spheroids to be of uniform 
density. User must set NC0MP=2 and provide dielectric function files for each spheroid. 

TWOELL - Two touching, homogeneous, isotropic ellipsoids, with distinct com- 
positions. 

SHPARl, SHPAR2, SHPAR3=x-length/c?, y-length/d, z-length/d of one ellipsoid. The two 
ellipsoids have identical shape, size, and orientation, but distinct dielectric functions. The 
line connecting ellipsoid centers is along the a;-axis in the TF. Target axes ai = (1,0,0) 
[along line connecting ellipsoids] and a2 = (0, 1, 0). User must set NC0MP=2 and provide 
dielectric function file names for both ellipsoids. Ellipsoids are in order of increasing x: 
first dielectric function is for ellipsoid with center at negative x, second dielectric function 
for ellipsoid with center at positive x. 

TWOAEL - Two touching, homogeneous, anisotropic ellipsoids, with distinct com- 
positions. 

Geometry as for TWOELL; SHPARl, SHPAR2, SHPAR3 have same meanings as for TWOELL. Tar- 
get axes ai = (1, 0, 0) and a2 = (0, 1, 0) in the TF. It is assumed that (for both ellipsoids) 
the dielectric tensor is diagonal in the TF. User must set NC0MP=6 and provide xx, yy, zz 
components of dielectric tensor for first ellipsoid, and xx, yy, zz components of dielectric 

tensor for second ellipsoid (ellipsoids are in order of increasing x). 

THRELL - Three touching homogeneous, isotropic ellipsoids of equal size and 
orientation, but distinct compositions. 

SHPARl, SHPAR2, SHPAR3 have same meaning as for TWOELL. Line connecting ellipsoid 
centers is parallel to x axis. Target axis ai = (1,0,0) along line of ellipsoid centers, 
a2 = (0, 1, 0). User must set NC0MP=3 and provide (isotropic) dielectric functions for first, 
second, and third ellipsoid. 

THRAEL - Three touching homogeneous, anisotropic ellipsoids with same size 
and orientation but distinct dielectric tensors. 

SHPARl, SHPAR2, SHPAR3 have same meanings as for THRELL. Target axis ai = (1,0,0) 

along line of ellipsoid centers, a2 = (0, 1,0). It is assumed that dielectric tensors are all 
diagonal in the TF. User must set NC0MP=9 and provide xx, yy, zz elements of dielectric 
tensor for first ellipsoid, xx, yy, zz elements for second ellipsoid, and xx, yy, zz elements 
for third ellipsoid (ellipsoids are in order of increasing x). 

BLOCKS - Homogeneous target constructed from cubic "blocks" . 

Number and location of blocks are specified in separate file blocks. par with following 

structure: 

one line of comments (may be blank) 

PRIN (= or 1 - see below) 
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N (= number of blocks) 

B (= width/d of one block) 

X y z {= position of 1st block in units of Bd) 

X y z {= position of 2nd block) 

X y z = position of Nth block 
If PRIN=0, then ai = (1, 0, 0), a2 = (0, 1,0). If PRIN=1, then ai and a2 are set to principal 
axes with largest and second largest moments of inertia, assuming target to be of uniform 
density. User must set NC0MP=1. 

• DW1996 - 13 block target used by Draine &: Weingartner (1996). 

Single, isotropic material. Target geometry was used in study by Draine & Weingartner 
(1996) of radiative torques on irregular grains, ai and a.2 are principal axes with largest 
and second-largest moments of inertia. User must set NC0MP=1. Target size is controlled 
by shape parameter SHPAR(l) = width of one block in lattice units. 

• NSPHER - Multisphere target consisting of the union of N spheres of single 
isotropic material. 

Spheres may overlap if desired. The relative locations and sizes of these spheres arc defined 
in an external file, whose name (enclosed in quotes) is passed through ddscat.par. The 
length of the file name should not exceed 13 characters. The pertinent line in ddscat .par 

should read 

SHPAR(l) SHPAR(2) '/ifename' (quotes must be used) 

where SHPAR(l) = target diameter in x direction (in Target Frame) in units of d 
SHPAR(2)= to have ai = (1,0,0), a2 = (0, 1,0) in Target Frame. 

SHPAR(2)= 1 to use principal axes of moment of inertia tensor for ai (largest /) and 02 
(intermediate I). 

filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined 
by parameter SHPAR(l), which is the extent of the multisphere target in the a;-direction, 
in units of the lattice spacing d. The file 'filename ' should have the following structure: 

A'^ (= number of spheres) 

one line of comments (may be blank) 

Xi yi z\ a\ (arb. units) 

Xi yi Z2 a2 (arb. units) 

xn yN zn a,N (arb. imits) 
where the coordinates of the center, and aj is the radius of sphere j. 

Note that xj, yj, zj, aj {j = 1,...,N) establish only the shape of the AT— sphere target. 
For instance, a target consisting of two touching spheres with the line between centers 
parallel to the x axis could equally well be described by lines 3 and 4 being 

0.5 

1 0.5 

or 

1 

2 1 

The actual size (in physical units) is set by the value of UeS specified in ddscat.par, 
where, as always, agg = (3y/47r)^/^, where V is the volume of material in the target. 
User must set NCDMP=1. 

• PRISMS - Triangulcir prism of homogeneous, isotropic material. 
SHPARl, SHPAR2, SHPAR3, SHPAR4 = a/d, b/a, c/a, L/a 
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The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the 
lattice spacing. The triangular cross-section has interior angles a, (3, 7 (opposite sides a, 6, 
c) given by cosa = {b'^ + — a^)/2bc, cos/? = (a^ + c^ — 6^)/2ac, COS7 = {a"^ + b'^ — c^)/2ab. 
In the Target Frame, the prism axis is in the ic direction, the normal to the rectangular 
face of width a is (0,1,0), the normal to the rectangular face of width b is (0, — COS7, sin 7), 
and the normal to the rectangular face of width c is (0, — cos 7, — sin 7). 

• FRMFIL - List of dipole locations and "compositions" obtained from a file. 

This option causes DDSCAT to read the target geometry information from a file shape . dat 
instead of automatically generating one of the geometries listed above. The shape . dat file 
is read by routine REASHP (file reashp.f). The user can customize REASHP as needed to 
conform to the manner in which the target geometry is stored in file shape . dat. However, 
as supplied, REASHP expects the file shape . dat to have the following structure: 

— one line containing a description; the first 67 characters will be read and printed in 
various output statements. 

— N = number of dipoles in target 

— aix aiy aiz — x,y,z components of ai 

— a2x a2y a2z = x,y,z components of a2 

— (line containing comments) 

- dummy IXYZ(1,1) IXYZ(1,2) IXYZ(1,3) ICDMP(1,1) IC0MP(1,2) IC0MP(1,3) 

- dummy IXYZ (2,1) IXYZ(2,2) IXYZ(2,3) ICDMP(2,1) IC0MP(2,2) IC0MP(2,3) 

- dummy IXYZ (3,1) IXYZ(3,2) IXYZ(3,3) IC0MP(3,1) IC0MP(3,2) IC0MP(3,3) 

- dummy ITild ,1) IXYZ(J,2) IXYZ(J,3) ICDMP(J,1) IC0MP(J,2) IC0MP(J,3) 

- dummy IXYZ (Yi ,1) IXYZ(N,2) IXYZ(N,3) ICDMP(N,1) IC0MP(N,2) IC0MP(N,3) 

where the number dummy is ignored, IXYZ(J,l-3) ~ x/d, y/d, z/d (TF) for J-th dipole, 
and IC0MP(J,l-3) = composition identifier to specify dielectric function appropriate for x,y,z 
directions at location of J-th dipole. For the common case of a single isotropic material, 
IC0MP(J,l-3) =111. 

The user should be able to easily modify these routines, or write new routines, to generate 
targets with other geometries. The user should first examine the routine target . f and modify 
it to call any new target generation routines desired. Alternatively, targets may be generated 
separately, and the target description (locations of dipoles and "composition" corresponding to 
x,y,z dielectric properties at each dipole site) read in from a file by invoking the option FRMFIL 
in ddscat . f . 

It is often desirable to be able to run the target generation routines without running the 
entire DDSCAT code. We have therefore provided a program CALLTARGET which allows the 
user to generate targets interactively; to create this executable just typef^ 
make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may 



need to edit the code to change the device number IDVOUT as for DDSCAT (see §3.1 above). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target. out is the same as 



Non-Unix sites: The source code for CALLTARGET is in the file MISC. FOR. You must compile and link this to 
ERRMSG, GETSET, LAPACKBLAS, LAPACKSUBS, PRINAXIS, REASHP, TAR2EL, TAR2SP, TAR3EL, TARBLOCKS, TARCEL, TARCYL, 
TARELL, TARGET, TARHEX, TARNSP, TARREC, TARTET, and WRIMSG. The source code for ERRMSG is in SRC2.F0R, GETSET is 
in SRC5 . FOR, and the rest are in SRC8 . FOR and SRC9 . FOR . 
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the format of the shape.dat files read if option FRMFIL is used (see above). Therefore you can 
simply 

mv target. out shape.dat 
and then use DDSCAT with the option FRMFIL in order to input a target shape generated by 
CALLTARGET. 

20 Scattering Directions 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix 
are reported in the output files vxxryyk.zzz.sc3i . The scattering direction is specified through 
angles 9s and (j)s (not to be confused with the angles © and $ which specify the orientation of 
the target relative to the incident radiation!). 

The scattering angle 9s is simply the angle between the incident beam (along direction x) 
and the scattered beam {9s = for forward scattering, dg = 180° for backscattering). 

The scattering angle specifies the orientation of the "scattering plane" relative to the x, y 
plane in the Lab Frame. When = the scattering plane is assumed to coincide with the x, y 
plane. When <ps = 90° the scattering plane is assumed to coincide with the x, z plane. Within 
the scattering plane the scattering directions are specified by < < 180°. 

Scattering directions for which the scattering properties are to be calculated are set in the 
file ddscat.par by specifying one or more scattering planes (determined by the value of (j)s) 
and for each scattering plane, the number and range of 9s values. The only limitation is that 
the number of scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code 
as distributed it is set to MXSCA=1000). 

21 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the 
X axis. DDSCAT. 5a allows the user to speciiy a general elliptical polarization for the incident 
radiation, by specifying the (complex) polarization vector eoi- The orthonormal polarization 
state eo2 = x x e^^ is generated automatically if ddscat .par specifies I0RTH=2. 

For incident linear polarization, one can simply set eoi = y by specifying (0,0) (1,0) 
(0,0) in ddscat.par; then eo2 = z For polarization mode eoi to correspond to right-handed 
circular polarization, set eoi = (y + i'z)/\/2 by specifying (0,0) (1,0) (0,1) in ddscat.par 
(DDSCAT. 5a automatically takes care of the normalization of eoi); then eo2 = («y + z)/^2, 
corresponding to left-handed circular polarization. 

22 Averaging over Scattering: g{l) = {cosOg), etc. 

DDSCAT automatically carries out numerical integration of various scattering properties. 
In particular, it calculates the mean value g{l) = {cos9s) for the scattered intensity for each 
incident polarization state. This is accomplished by evaluating the scattered intensity for ICTHM 
different values of Og (including 6s =0 and 6s = tt), and taking a weighted sum. For each value 
of 9s except and tt, the scattering intensity is evaluated for IPHM different values of the 
scattering angle . The angular integration over 9s is accomplished using Simpson's rule (with 
uniform intervals in cos^g), so ICTHM should be an odd number. The angular integration over 
4>s is accomplished by sampling uniformly in with uniform weighting; IPHM can be either 
even or odd. 

The following quantities are evaluated by this angular integration: 
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• g = (cos^s)^ + (sin 0s cos0s)y + (sinflg sin0s)z (see j|T^); 

• Qr (see §|I|). 

It is important that the user recognize that accurate evaluation of these angular averages requires 
adequate sampling over scattering angles. For small values of the size parameter x — 27racff/A, 
the angular distribution of scattered radiation has a dipolar character and the sampling in 9^ and 
(ps does not need to be very fine, so ICTHM and IPHM need not be large. For larger values of the 
size parameter x, however, higher multipoles in the scattered radiation field become important, 
and finer sampling in 9g and <f>g is required. We do not have any foolproof prescription to offer, 
since the scattering pattern will depend upon the target geometry and dielectric constant in 
addition to overall size parameter. However, as a very rough guide, we suggest that the user 
specify values of ICTHM and IPHM satisfying 

ICTHM > 5(1 + 2;) , (30) 
IPHM > 2(1 + 2;) . (31) 

The sample ddscat.par file supplied has ICTHM — 33 and IPHM = 12; the above criteria would 
suggest that this would be suitable for 2; < 5. 

The cpu time required for evaluation of these angular averages is proportional to [2 + 
IPHM(ICTHM — 2)]. Since the computational time spent in evaluating these angular integrals 
can be a significant part of the total, it is important to choose values of ICTHM and IPHM which 
will provide a suitable balance between accuracy (in this part of the overall calculation) and 
cpu time. 

Within one scattering plane, the scattered intensity tends to have approximately (1 + x) 
peaks for < 6*3 < tt, so that the above prescription for ICTHM would have at least 5 sampling 
points per maximum. The angular distribution over (ps is usually not as structured as that over 
6s so we suggest that IPHM need not be as large as ICTHM. We have refrained from "hard-wiring" 
the values of ICTHM and IPHM because we are not confident of the reliability of the recommended 
criteria ( |3C| , |3l| ) - it is up to the user to specify appropriate values of ICTHM and IPHM according 
to the requirements of the problem being addressed. 



23 Mueller Matrix for Scattering in Selected Directions 

23.1 Two Orthogonal Incident Polarizations (I0RTH=2) 

DDSCAT.5a internally computes the scattering properties of the dipole array in terms of a 
complex scattering matrix fmi{0s,4's) (Draine 1988), where index / = 1,2 denotes the incident 
polarization states, m = 1,2 denotes the scattered polarization state, and 9s,(f>s specify the 
scattering direction. Normally DDSCAT is used with I0RTH=2 in ddscat.par, so that the 
scattering problem will be solved for both incident polarization states (^ = 1 and 2); in this 
subsection it will be assumed that this is the case. 

Incident polarization states I = 1,2 correspond to polarization states egi, 692; recall that 
polarization state eoi is user-specified, and §02 = x x e^^. Scattered polarization state m — 1 
corresponds to linear polarization of the scattered wave parallel to the scattering plane (ei = 

= 0s) and m = 2 corresponds to linear polarization perpendicular to the scattering plane (in 
the +0S direction). The scattering matrix fmi was defined (Draine 1988) so that the scattered 
electric field Eg is related to the incident electric field Ei(0) at the origin (where the target is 
assumed to be located) by 

f l]s-0s \_ exp(^k ■ r) / /n /12 \ / E,(0) • e*o, \ , . 

kr \hi /22 A E.(0)-e52 ; • ^''^^ 



33 



The 2x2 complex amplitude scattering matrix (with elements ^i, 52, 5*3, and ^4) is defined so 
that (see Bohren & Huffman 1983) 



E, -Os \ _ cxp(ik • r) / ^2 S3\ f E,(0) • e 



• 0s / -ikr \ Si Si J \ Ei(0)-e, 



(33) 



where ej||, e^j^ are (real) unit vectors for incident polarization parallel and perpendicular to the 
scattering plane (with the customary definition of e^^ = e^n x x). 
From (|3^ , |33| ) we may write 



^2 S,\f E,(0) -6,11 \_^,f hi fi2 \ ( E,(0) -e*! 
^4 Si \ E,(0) •e.i i I -/21 -/22 j I E,(0) -eS^ 



(34) 



Let 



(35) 
(36) 
(37) 
(38) 

Note that since eoi,eo2 could be complex (i.e., elliptical polarization), the quantities a,b,c,d 
are complex. Then 

a b \ f y 

c d 



a — 


^01 


•y 


b = 


^01 


• z 


c = 


^02 


•y 


d = 


^02 


• z 



(39) 



and eq. (|3j) can be written 

^2 S3\f E,(0) • e,|| \ / -/n -/12 \{ a b\f E,(0) • y 
Si Si ){E,{0)-e,^ ) /21 /22 )\c d ){E,{0)-i 

The incident polarization states ej|| and e^j^ are related to y, z by 

y \ ^ ( cos 0s sin 0s \ f e^ii 
z I \ sin 03 -cos (As / V ei± 



(40) 



(41) 



substituting (^TJ) into ( pO| ) we obtain 

^2 53 \ / E,(0) • e,|| \ / -fii -fi2 \f a 6 \ / cos 0s sin 0s \ f E,(0) • y 
^4 Si J\ E,(0) -e,^ ) \ f2i f22 )\c d )\ sin0s -cos0s ) \ E,(0) • z 

(42) 

Eq. ( ^ ) must be true for all Ei(0), so we obtain an expression for the complex scattering 
amplitude matrix in terms of the fmi '■ 



S2 '^3 \ ^-f -III -/i2 \ ( ^ 6 \ / cos 0s sin0. 
Si Si ) V /21 f22 J \ c d ) \ sin 0s -cost 



(43) 



This provides the 4 equations used in subroutine GETMUELLER to compute the amplitude scat- 
tering matrix elements: 

51 = — i [/2i(6cos0s — asin0s) + /22(rfcos0s — ccos0s)] , (44) 

52 = [/ii(acos0s + 6sin0s) + /i2(ccos0s + dsin0s)] , (45) 

53 = i [/ii(5cos0s - asin0s) + /i2(dcos0s - csin0s)] , (46) 
'S'4 = i [/2i(acos0s + 6sin0s) + /22(ccos0s + (isin0s)] . (47) 
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It is both convenient and customary to characterize both incident and scattered radiation by 
4 "Stokes parameters" - the elements of the "Stokes vector" . There are diflFerent conventions 
in the hterature; we adhere to the definitions of the Stokes vector {I,Q,U,V) adopted in the 
excellent treatise by Bohren & Huffman (1983), to which the reader is referred for further detail. 
Here are some examples of Stokes vectors (/, Q,U,V): 

(1,0,0,0)/ : impolarized light. 

(1, 1, 0, 0)/ : 100% linearly polarized with E parallel to the scattering plane; 

(1, —1, 0, 0)7 : 100% linearly polarized with E perpendicular to the scattering plane; 

(1, 0, 1, 0)7 : 100% linearly polarized with E at +45° relative to the scattering plane; 

(1, 0, —1, 0)7 : 100% linearly polarized with E at -45° relative to the scattering plane; 

(1,0,0, 1)7 : 100% right circular polarization {i.e., negative helicity); 

(1,0,0,-1)7 : 100% left circular polarization {i.e., positive helicity). 

It is convenient to describe the scattering properties in terms of the 4x4 Mueller matrix relating 
the Stokes parameters {Ii,Qi, Ui, Vi) and (7s, Qs, Us, Vs) of the incident and scattered radiation: 







f Sn 


S'12 


Sl3 


Sl4 \ 




( \ 




1 


S21 


S22 


S23 


S24 




Qi 


Us 




S31 


S32 


S33 


S34 




u. 






\ ^41 


842 


843 


S44 1 







(48) 



This 4x4 Mueller matrix for scattering by a single particle is also referred to as the "scattering 
matrix" , and sometimes as the "phase matrix" . 

Once the amplitude scattering matrix elements Si, S2, S3, and S4 have been obtained, the 
Mueller matrix elements can be computed (Bohren & Huffman 1983): 



Sii 




{\Sif + \S2f + \S3f 


+ \S4\')/2 


S12 




{\S2\^-\Sl\' + \S4\^ 


-\S3\')/2 


Sl3 




Re (S'25'3 + 5'iS'4 ) , 




Sl4 




lm{S2S3-SiSl) , 




S2I 




{\S2f-\Sif + \S3f 




S22 




{\Slf + \S2\^-\S3f 


-\S4\')/2 


S23 




Re{S2S^- SiSl) , 




S24 




Im {S2S^ + SiSl) , 




S3I 




Re{S2Sl + SiS;) , 




S32 




Re{S2Sl-SiS;) , 




S33 




Re{SiS; + S3Sl) , 




S34 




Im{S2Sl + S4S^) , 




S41 




lm{S4S^+SiS;) , 




S42 




Im {S4S2 — SiS^) , 




S43 




Im {S1S2 - S3SI) , 




S44 




Re (S'iS'2 — >S'35'4) 





(49) 



These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which 
handles output of scattering properties. As delivered, WRITESCA writes out 6 selected elements: 
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Sii, S21, S31, S41 (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), 6*12, and 5*13. In addition, WRITESCA writes out the linear 
polarization P of the scattered light for incident unpolarized light: 

Oil 

Of course, other elements Sij may be of interest. It is relatively straightforward for the user 
to modify subroutine WRITESCA to write out whatever elements of the Mueller matrix (or the 
scattering amplitude matrix) are desired. 



23.2 One Incident Polarization State Only (I0RTH=1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state 
- for example, when each solution is very time-consuming, and the target is known to have 
some symmetry so that solving for a single incident polarization state may be sufficient for the 
required purpose. In this case, set IDRTH=1 in ddscat .par. 

When I0RTH=1, only /n and /21 are available; hence, DDSCAT cannot automatically 
generate the Mueller matrix elements. In this case, the output routine WRITESCA writes out the 
quantities |/iip, I/21P, ^^{fiifii)^ Iin(/ii/2i) foi' ^&ch of the scattering directions. 

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 
even if ddscat. par specified I0RTH=1: this is because when more than one value of the target 
orientation angle $ is required, there is no additional "cost" to solve the scattering problem 
for the second incident polarization state, since when solutions are available for two orthogonal 
states for some particular target orientation, the solution may be obtained for another target 
orientation differing only in the value of $ by appropriate linear combinations of these solutions. 
Hence we may as well solve the "complete" scattering problem so that we can compute the 
complete Mueller matrix. 



24 Graphics and Postprocessing 
24.1 IDL 

At present, we do not offer a comprehensive package for DDSCAT data postprocessing and 
graphical display in IDL. However, there are several developments worth mentioning: First, 



we offer several output capabilities from within DDSCAT: ASCII (sec § |l0.lD , FORTRAN 
unformatted binary (see § 10.2 ), and netCDF portable binary (see § |10.3| ). Second, we offer 



several skeleton IDL utilities: 

• bhmie . pro is our translation to IDL of the popular Bohren-Huffman code which calculates 
efficiencies for spherical particles using Mie theory. 

• readbin . pro reads the FORTRAN unformatted binary file written by routine writebin . f . 
The variables are stored in a common block. 

• readnet .pro reads NetCDF portable binary file and should be the method of choice for 
IDL users. It offers random data access. 

• mie. pro is an example of an interface to binary files and the Bohrcn-Huffman code. It 
plots a comparison of DDSCAT results with scattering by equivalent radius spheres. 

At present the IDL code is experimental. 
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25 Miscellanea 



Additional source code, refractive index files, etc., contributed by users will be located in the 
directory DDA/misc. These routines and files should be considered to be not supported by 
Draine and Flatau - caveat receptor! These routines and files should be accompanied by enough 
information (e.g., comments in source code) to explain their use. 



26 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure 
of the ddscat .par file is intended to be simple and suggestive so that, after reading the above 
notes once, the user may not have to refer to them again. 

The file rel_notes in DDA/doc lists known bugs in DDSCAT. 5a . An up-to-date version 
will be maintained in the directory where the publicly available source code is located (see the 
README file). Users are encouraged to provide B. T. Draine (draine@astro.princeton.edu) 
with their email address; bug reports, and any new releases of DDSCAT, will be made known 
to those who do! 

P. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" 



with URL http://atol.ucsd.edu/--pflatau. The SCATTERLIB Internet site is a library of 



light scattering codes. Emphasis is on providing source codes (mostly FORTRAN). However, 
other information related to scattering on spherical and non-spherical particles is collected: an 
extensive list of references to light scattering methods, refractive index, etc. This URL page 
contains a section on the discrete dipole approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. If you 
wish to cite this User Guide, we suggest the following citation: 

Draine, B.T., & Flatau, P. J. 2000, "User Guide for the Discrete Dipole Approximation Code 
DDSCAT (Version 5alO)" , |http://xxx.lanl.gov/abs/astro-ph/0008151v2| 



Finally, the authors have one special request: We would very much appreciate preprints and 
(especially!) reprints of any papers which make use of DDSCAT! 
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A Understanding and Modifying ddscat.par 



In order to use DDSCAT to perform the specific calculations of interest to you, it will be 
necessary to modify the ddscat .par file. Here wc list the sample ddscat .par file, followed by a 
discussion of how to modify this file as needed. Note that all numerical input data in DDSCAT 
is read with free-format READ(IDEV,*) . . . statements. Therefore you do not need to worry 
about the precise format in which integer or floating point numbers arc entered on a line. The 
crucial thing is that lines in ddscat .par containing numerical data have the correct number of 
data entries, with any informational comments appearing after the numerical data on a given 
line. 

' =================== Parameter file ===================' 

PRELIMINARIES ****' 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGST'= CMDS0L*6 (PBCGST, PETRKP) — select solution method 
'GPFAFT'= CMETHD*6 (GPFAFT, BRENNR, TMPRTN, CONVEX) 
'LATTDR'= CALPHA*6 (LATTDR, LDRISO, G0BR88, DRAI88) 
'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 
'NOTCDF'= CNETFLAG (ALLCDF, ORICDF, NOTCDF) 

' RCTNGL ' = CSHAPE*6 (FRMFIL , ELLIPS , CYLNDR , RCTNGL , HEXGON , TETRAH , UNICYL .UNIELL) 

8 6 4= shape parameters PARI, PAR2, PARS 

1 = NCOMP = number of dielectric materials 

'TABLES'= CDIEL*6 (TABLES, H20ICE,H20LIQ; if TABLES, then filenames follow...) 
'diel.tab' 

'**** CONJUGATE GRADIENT DEFINITIONS ****' 

= INIT (TO BEGIN WITH |X0> = 0) 

l.OOe-5 = ERR = MAX ALLOWED (NORM OF |G>=AC|E>-ACA|X>)/(NORM OF AC|E>) 
'**** ANGLES FOR CALCULATION OF CSCA, G' 

33 = ICTHM (number of theta values for evaluation of Csca and g) 

12 = IPHM (number of phi values for evaluation of Csca and g) 

'**** Wavelengths (micron) ****' 

6.283185 6.283185 1 'INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

1. 1. 1 'LIN' = eff. radii (first, last, how many, how=LIN,INV,LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x Eixis) 

2 = lORTH (=1 to do only pol. state eOl; =2 to also do orth. pol. state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each teirget orient. 
'*+** Prescribe Target Rotations ***♦' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify Scattered Directions ***** 

0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180 . 30 = phi , ... for plane B 
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Lines Comments 

1-2 comment lines - no need to change. 

3 NOTORQ if torque calculation is not required; 
DOTORQ if torque calculation is required. 

4 PBCGST is recommended; other option is PETRKP (see §|l^. 

5 GPFAFT recommended; other options are BRENNR, TMPRTN, CONVEX (§|l|). 

6 LATTDR (LATTice Dispersion Relation) is rec ommended (§|ll|). 

7 ALLBIN for unformatted binary dump (§10.2; 

ORIBIN for unformatted binary dump of orientational averages only; 
NOTBIN for no unformatted binary output. 

8 ALLCDF for output in netCDF format (must have netCDF option enabled; cf. § 10. 3| ); 
ORICDF for orientational averages in netCDF format (must have netCDF option enabled). 
NOTCDF for no output in netCDF format; 

9 specify choice of target shape (see §|l^for description of options RCTNGL, ELLIPS, TETRAH, 

10 shape parameters SHPARl, SHPAR2, SHPAR3, ... (sec §|l|). 

11 number of different dielectric constant tables (§p^). 

12 name(s) of dielectric constant table(s) (one per line). 

13 comment line - no need to change. 

14 is recommended value of parameter INIT. 

15 ERR = error tolerance h: maximum allowed value of E — A'^ AP\/\A^ E\ [see eq.(|l8|)]. 

16 comment line - no need to change. 

17 ICTHM - number of Og values for angular averages (§p2[). 

18 IPHM - number of (j)s values for angular averages. 

19 comment line - no need to change. 

20 A - first, last, how many, how chosen. 

21 comment line - no need to change. 

22 aeff - first, last, how many, how chosen. 

23 comment line - no need to change. 

24 specify x,y,z components of (complex) incident polarization Gqi (§|2l|) 

25 lORTH = 1 to do one polarization state only; 

2 to do second (orthogonal) incident polarization as well. 

26 IWRKSC = to suppress writing of ".sea" files; 
2 to enable writing of ".sea" files. 

27 comment line - no need to change. 

28 /3 (see § p7| ) - first, last, how many . 

29 Q - first, last, how many. 

30 $ - first, last, how many. 

31 comment line - no need to change. 

32 0s for first scattering plane, 9s,min, Os,max, how many 9s values; 
33,... (j)s for 2nd,... scattering plane, ... 
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B yixxryyorl . avg Files 



The file wOOrOOori . avg contains the results for the first wavelength (wOO) and first target radius 
(rOO) averaged over orientations (or i. avg). The wOOrOOori . avg file generated by the sample 
calculation should look like the following: 

DDSCAT DDSCAT.5a.10 [00.11.03] 

TARGET Rectangular prism; NX,NY,NZ= 8 6 4 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

NATO = 192 = number of dipoles 
AEFF= 1.00000 = effective radius (physical units) 
WAVE= 6.28319 = wavelength (physical units) 

K*A= 1.00000 = 2*pi*aeff /lambda 
n= ( 1.3300 , 0.0100), epsilon= ( 1.7688 , 0.0266) for material 1 

TOL= l.OOOE-05 = error tolerance for CCG method 
ICTHM= 33 = theta values used in comp. of Qsca.g 
IPHIM= 12 = phi values used in comp. of Qsca.g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.27942 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, 0. 00000) =inc.pol.vec. 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, 0. 00000) =inc.pol.vec. 

0.000 0.000 = beta_min, beta_max ; NBETA = 1 

0.000 90.000 = theta_min, theta_max; NTHETA= 3 

0.000 0.000 = phi_min, phi_max ; NPHI = 1 
Results averaged over 3 target orientations 

and 2 incident polarizations 



in LF 
in LF 



J0=1 
J0=2 
mean 
Qpol^ 

J0=1 
J0=2 
mean 



Qext Qabs Qsca 

1.330E-01 3.269E-02 1.003E-01 
9.063E-02 2.397E-02 6.666E-02 
1.118E-01 2.833E-02 8.347E-02 
4.234E-02 
Qsca*g(l) 
2.352E-02 
1.781E-02 
2.066E-02 



Qsca*g(2) 
1.126E-03 
2 . 768E-03 
1 . 947E-03 



g=<cos> Qbk Qpha 
2.345E-01 5.552E-03 4.836E-01 
2.673E-01 3.430E-03 4.159E-01 
2.476E-01 4.491E-03 4.497E-01 
dQpha= 6.768E-02 



Qsca*g(3) 
-7. 190E-10 
3.244E-11 
-3.433E-10 



** Mueller matrix elements for selected scattering directions ** 



theta 


phi 




S_ll 






S_21 






S_31 






S_41 






S_12 






S_13 




Pol. 














5 


167E- 


-02 


7 


916E- 


-03 


9 


582E- 


■11 


1 


267E- 


■11 


7 


916E- 


-03 


9 


239E- 


11 





15320 


30 











4 


046E- 


-02 


4 


518E- 


-04 


2 


540E- 


-11 


4 


619E- 


■11 


4 


518E- 


-04 


2 


339E- 


11 





01117 


60 











2 


169E- 


-02 


-1 


023E- 


-02 


8 


752E- 


-11 


5 


628E- 


■11 


-1 


023E- 


-02 


6 


967E- 


■11 





47184 


90 











1 


193E- 


-02 


-1 


192E- 


-02 


2 


068E- 


■11 


6 


778E- 


■11 


-1 


192E- 


-02 


-1 


481E- 


11 





99888 


120 











1 


170E- 


-02 


-5 


702E- 


-03 


-4 


564E- 


■11 


3 


005E- 


■11 


-5 


702E- 


-03 


5 


392E- 


11 





48743 


150 











1 


403E- 


-02 


9 


758E- 


-04 


-2 


008E- 


■11 


-1 


155E- 


11 


9 


758E- 


-04 


1 


665E- 


11 





06956 


180 











1 


411E- 


-02 


3 


333E- 


-03 


-1 


390E- 


■11 


1 


454E- 


■11 


3 


333E- 


-03 


1 


643E- 


11 





23621 








90 





5 


167E- 


-02 


-7 


916E- 


-03 


-7 


574E- 


■10 


4 


831E- 


■11 


-7 


916E- 


-03 


-7 


428E- 


10 





15320 


30 





90 





4 


487E- 


-02 


-1 


280E- 


-02 


-5 


036E- 


■04 


8 


979E- 


■05 


-1 


281E- 


-02 


-3 


843E- 


04 





28550 


60 





90 





3 


040E- 


-02 


-2 


062E- 


-02 


-7 


990E- 


■04 


1 


643E- 


■04 


-2 


063E- 


-02 


-3 


773E- 


04 





67857 


90 





90 





1 


991E- 


-02 


-1 


989E- 


-02 


-7 


361E- 


■04 


1 


827E- 


■04 


-1 


991E- 


-02 


-3 


807E- 


05 





99926 


120 





90 





1 


626E- 


-02 


-1 


175E- 


-02 


-4 


415E- 


■04 


1 


359E- 


■04 


-1 


176E- 


-02 


1 


634E- 


■04 





72285 


150 





90 





1 


478E- 


-02 


-5 


293E- 


-03 


-1 


716E- 


■04 


6 


554E- 


■05 


-5 


295E- 


-03 


1 


131E- 


■04 





35833 


180 





90 





1 


411E- 


-02 


-3 


333E- 


-03 


3 


067E- 


■10 


-2 


982E- 


■11 


-3 


333E- 


-03 


-3 


009E- 


10 





23621 
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C Mixxryykzzz. sea Files 



The wOOrOOkOOO. sea file contains the results for the first wavelength (wOO), first target radius 
(rOO), and first orientation (kOOO). The wOOrOOkOOO . sea file created by the sample calculation 
should look like the following: 

DDSCAT DDSCAT.5a.10 [00.11.03] 

TARGET Rectangular prism; NX,NY,NZ= 8 6 4 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

NATO = 192 = number of dipoles 
AEFF= 1.00000 = effective radius (physical units) 
WAVE= 6.28319 = wavelength (physical units) 

K*A= 1.00000 = 2*pi*aeff /lambda 
n= ( 1.3300 , 0.0100), epsilon= ( 1.7688 , 0.0266) for material 1 

TOL= l.OOOE-05 = error tolerance for CCG method 
ICTHM= 33 = theta values used in comp. of Qsca.g 
IPHIM= 12 = phi values used in comp. of Qsca.g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.27942 0.00000 0.00000) = k vector (latt. units) in TF 

( 0.00000, 0.00000)( 1.00000, 0.00000)( 0.00000, 0. 00000) =inc.pol.vec. 1 in TF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, 0. 00000) =inc.pol.vec. 2 in TF 
BETA = 0.000 = rotation of target around Al 
THETA= 0.000 = angle between Al and k 
PHI = 0.000 = rotation of Al aromid k 





Qext 


Qabs 


Qsca g=<cos> 




Qbk 




Qpha 




J0=1 


l.llOE-01 


3.028E-02 8 


077E-02 3.504E-01 


1 


858E-03 


4 


668E- 


-01 


J0=2 


8.651E-02 


2.441E-02 6 


209E-02 3.650E-01 


1 


362E-03 


4 


197E- 


-01 


mean 


9.878E-02 


2.735E-02 7 


143E-02 3.567E-01 


1 


610E-03 


4 


432E- 


-01 


Qpol= 


= 2.454E-02 
Qsca*g(l) 


Qsca*g(2) 


Qsca*g(3) 




dqpha= 


4 


711E- 


-02 


JQ=1 


2.830E-02 


9.829E-10 - 


-2.070E-09 












J0=2 


2 . 266E-02 


1 . 827E-09 


9.086E-11 












mean 


2 . 548E-02 


1.405E-09 - 


-9.896E-10 













** Mueller matrix elements for selected scattering directions ** 



theta 


phi 




S_ll 






S_21 






S_31 






S_41 






S_12 






S_13 




Pol. 














4 


987E- 


-02 


5 


372E- 


-03 


8 


485E- 


11 


-1 


807E- 


11 


5 


372E- 


-03 


8 


482E- 


11 





10772 


30 











4 


040E- 


-02 


-9 


275E- 


-04 


-4 


151E- 


11 


1 


973E- 


11 


-9 


275E- 


-04 


-4 


057E- 


11 





02296 


60 











2 


191E- 


-02 


-1 


060E- 


-02 


1 


292E- 


10 


1 


107E- 


11 


-1 


060E- 


-02 


1 


205E- 


10 





48383 


90 











1 


058E- 


-02 


-1 


056E- 


-02 


-2 


798E- 


11 


7 


613E- 


11 


-1 


056E- 


-02 


-1 


444E- 


10 





99794 


120 











7 


506E- 


-03 


-3 


995E- 


-03 


3 


887E- 


11 


-3 


829E- 


11 


-3 


995E- 


-03 


-4 


776E- 


12 





53222 


150 











5 


921E- 


-03 


-9 


862E- 


-06 


1 


709E- 


11 


-8 


195E- 


13 


-9 


862E- 


-06 


-1 


717E- 


11 





00167 


180 











5 


058E- 


-03 


7 


791E- 


-04 


6 


177E- 


12 


-8 


607E- 


12 


7 


791E- 


-04 


-4 


347E- 


12 





15403 








90 





4 


987E- 


-02 


-5 


372E- 


-03 


-4 


752E- 


10 


4 


820E- 


11 


-5 


372E- 


-03 


-5 


005E- 


10 





10772 


30 





90 





4 


285E- 


-02 


-1 


032E- 


-02 


-5 


283E- 


10 


-1 


079E- 


11 


-1 


032E- 


-02 


-5 


098E- 


10 





24082 


60 





90 





2 


736E- 


-02 


-1 


773E- 


-02 


-2 


169E- 


10 


7 


855E- 


11 


-1 


773E- 


-02 


-2 


703E- 


10 





64798 


90 





90 





1 


540E- 


-02 


-1 


539E- 


-02 


-2 


506E- 


11 


2 


845E- 


11 


-1 


539E- 


-02 


-1 


588E- 


10 





99942 


120 





90 





9 


786E- 


-03 


-6 


745E- 


-03 


1 


827E- 


11 


3 


636E- 


12 


-6 


745E- 


-03 


-7 


262E- 


11 





68918 


150 





90 





6 


389E- 


-03 


-1 


820E- 


-03 


4 


853E- 


11 


2 


225E- 


11 


-1 


820E- 


-03 


-5 


195E- 


11 





28491 


180 





90 





5 


058E- 


-03 


-7 


79 lE- 


-04 


5 


552E- 


11 


2 


060E- 


12 


-7 


791E- 


-04 


-5 


362E- 


11 





15403 
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